EGR 511



EGR 599 Advanced Engineering Math II _____________________

LAST NAME, FIRST

Problem set #6

1. (P. Chapra 17.6) Use least-squares regression to fit a straight line, y = a + bx, to

|x |2 |3 |4 |7 |8 |9 |5 |5 |

|y |9 |6 |5 |10 |9 |11 |2 |3 |

(a) Along with the slope and intercept, compute the standard error of the estimate and the correlation coefficient. Plot the data and the straight line. Assess the fit.

(b) Recompute (a), but use polynomial to fit a parabola, y = a + bx + cx2, to the data. Compare the results with those of (a).

Solution

% Problem 5.9

x=[2 3 4 7 8 9 5 5];

y=[9 6 5 10 9 11 2 3];

n=length(x);

coef=polyfit(x,y,1);

ycal=polyval(coef,x);

S=sum((ycal-y).^2);

yave=mean(y);

Sdev=sum((y-yave).^2);

stde=sqrt(S/(n-2));

cor=sqrt(1-S/Sdev);

coef1=coef;

fprintf('y = a+bx, a = %8.5f, b = %8.5f\n',coef(2),coef(1))

fprintf('Standard error = %8.4f\n',stde)

fprintf('Correlation coefficient = %8.4f\n',cor)

coef=polyfit(x,y,2);

ycal=polyval(coef,x);

S=sum((ycal-y).^2);

yave=mean(y);

Sdev=sum((y-yave).^2);

stde=sqrt(S/(n-2));

cor=sqrt(1-S/Sdev);

fprintf('y = a+bx+cx^2, a = %8.5f, b = %8.5f, c = %8.5f\n',coef(3),coef(2),coef(1))

fprintf('Standard error = %8.4f\n',stde)

fprintf('Correlation coefficient = %8.4f\n',cor)

xmax=max(x);xmin=min(x);

x1=[xmin xmax];

y1=polyval(coef1,x1);

dx=(xmax-xmin)/50;

x2=xmin:dx:xmax;

y2=polyval(coef,x2);

plot(x,y,'o',x1,y1,x2,y2)

xlabel('x');ylabel('y')

>> s5d9

y = a+bx, a = 3.48955, b = 0.62985

Standard error = 3.2214

Correlation coefficient = 0.4589

y = a+bx+cx^2, a = 16.02696, b = -4.80692, c = 0.48894

Standard error = 1.9254

Correlation coefficient = 0.8474

[pic]

2. (P. Chapra 17.13) Given the data

x |5 |10 |15 |20 |25 |30 |35 |40 |45 |50 | |y |16 |25 |32 |33 |38 |36 |39 |40 |42 |42 | |

use least-squares regression to fit

(a) a straight line, y = a + bx,

(b) a power equation, y = axb,

(c) a saturation-growth-rate equation, y = a[pic], and

(d) a parabola, y = a + bx + cx2.

Plot the data along with all the curves in (a), (b), (c) and (d). Is any one of the curves superior? If so, justify.

(e) Use nonlinear regression to fit a saturation-growth-rate equation, y = a[pic], to the data.

Solution

% Problem 6.2

x=[5 10 15 20 25 30 35 40 45 50];

y=[16 25 32 33 38 36 39 40 42 42];

n=length(x);

coef=polyfit(x,y,1);

ycal=polyval(coef,x);

S=sum((ycal-y).^2);

yave=mean(y);

Sdev=sum((y-yave).^2);

stde=sqrt(S/(n-2));

cor=sqrt(1-S/Sdev);

coef1=coef;

fprintf('y = a+bx, a = %8.5f, b = %8.5f\n',coef(2),coef(1))

fprintf('Standard error = %8.4f\n',stde)

fprintf('Correlation coefficient = %8.4f\n',cor)

disp(' ')

xlog=log(x);ylog=log(y);

coef=polyfit(xlog,ylog,1);

a=exp(coef(2));b=coef(1);

ycal=a*x.^b;

S=sum((ycal-y).^2);

stde=sqrt(S/(n-2));

cor=sqrt(1-S/Sdev);

a2=a;b2=b;

fprintf('y = ax^b, a = %8.5f, b = %8.5f\n',a,b)

fprintf('Standard error = %8.4f\n',stde)

fprintf('Correlation coefficient = %8.4f\n',cor)

disp(' ')

xi=1.0./x;yi=1.0./y;

coef=polyfit(xi,yi,1);

a=1/coef(2);b=a*coef(1);

ycal=a*x./(b+x);

S=sum((ycal-y).^2);

stde=sqrt(S/(n-2));

cor=sqrt(1-S/Sdev);

a3=a;b3=b;

fprintf('y = ax/(b+x), a = %8.5f, b = %8.5f\n',a,b)

fprintf('Standard error = %8.4f\n',stde)

fprintf('Correlation coefficient = %8.4f\n',cor)

disp(' ')

coef=polyfit(x,y,2);

ycal=polyval(coef,x);

S=sum((ycal-y).^2);

stde=sqrt(S/(n-2));

cor=sqrt(1-S/Sdev);

fprintf('y = a+bx+cx^2, a = %8.5f, b = %8.5f, c = %8.5f\n',coef(3),coef(2),coef(1))

fprintf('Standard error = %8.4f\n',stde)

fprintf('Correlation coefficient = %8.4f\n',cor)

xmax=max(x);xmin=min(x);

x1=[xmin xmax];

y1=polyval(coef1,x1);

dx=(xmax-xmin)/50;

x2=xmin:dx:xmax;

y2=a2*x2.^b2;

y3=a3*x2./(b3+x2);

y4=polyval(coef,x2);

plot(x,y,'o',x1,y1,'k',x2,y2,'r:',x2,y3,'g-.',x2,y4,'b--')

xlabel('x');ylabel('y')

legend('data','y=a+bx','y=ax^b','y=ax/(b+x)','y=a+bx+cx^2',4)

>> s5d10

y = a+bx, a = 20.66667, b = 0.49576

Standard error = 3.7281

Correlation coefficient = 0.9056

y = ax^b, a = 9.65828, b = 0.39503

Standard error = 2.5472

Correlation coefficient = 0.9571

y = ax/(b+x), a = 52.24201, b = 11.14680

Standard error = 1.3251

Correlation coefficient = 0.9886

y = a+bx+cx^2, a = 12.16667, b = 1.34576, c = -0.01545

Standard error = 2.0115

Correlation coefficient = 0.9735

>>

(c) is best because standard error is smallest.

[pic]

(e) Use nonlinear regression to fit a saturation-growth-rate equation, y = a[pic], to the data.

% Gauss-Newton method

%

x=[5 10 15 20 25 30 35 40 45 50]';

y=[16 25 32 33 38 36 39 40 42 42]';

n=length(x);

a=50;b=10;

ymodel='a*x./(b+x)';

dyda='x./(b+x)';dydb='-a*x./(b+x).^2';

for i=1:2

Zj=[eval(dyda) eval(dydb)];

ZjTZj=Zj'*Zj

D=y-eval(ymodel)

ZjTD=Zj'*D

dA=ZjTZj\ZjTD

a=a+dA(1);

b=b+dA(2);

fprintf('Iteration #%g: a = %8.4f, b = %8.4f\n',i,a,b)

end

>> s6p2e

ZjTZj =

4.8471 -6.3875

-6.3875 9.8141

D =

-0.6667

0

2.0000

-0.3333

2.2857

-1.5000

0.1111

0

1.0909

0.3333

ZjTD =

2.5200

-3.3539

dA =

0.4887

-0.0237

Iteration #1: a = 50.4887, b = 9.9763

ZjTZj =

4.8531 -6.4641

-6.4641 10.0416

D =

-0.8562

-0.2743

1.6781

-0.6857

1.9122

-1.8889

-0.2897

-0.4101

0.6733

-0.0905

ZjTD =

-0.0016

0.0034

dA =

1.0e-003 *

0.8300

0.8694

Iteration #2: a = 50.4895, b = 9.9772

3. (Rowley, R.L.2, P. 7.2) Use an MC (Monte Carlo) method to compute the area of a regular hexagon inscribed in a circle of radius 5.

Solution

[pic]

We can formulate the problem in terms of the diagram shown above. By symmetry, we need only look at one quadrant. We will choose random coordinates (x, y) such that the distance from the origin is less than r = 5. The condition for “hits” is then when x ( r/2, y must be less than (r/2)(31/2; and when r/2 < x ( r, then y must be less than (r ( x)tan (/3 = (r ( x)(31/2. We can then divide the number of hits by the total number of trials and multiply by the area of the circle, (r2.

% Set 6, problem 3

trials=10000;

r = 5;

r2 = r/2; s3=sqrt(3); h = r2*sqrt(3);

% Area of right triangle with side r2 and h

% Area of hexagon is 12 times area of right triangle

Actual_area = 6*r2*h;

hits = 0;

for i=1:trials

x = r; y = r;

while sqrt(x*x+y*y)>r

x = r*rand(1);y = r*rand(1);

end

if x< r2

if y s6p3

trials = 10000 MC area = 65.1252 Actual area = 64.9519

4. (Rowley, R.L., P. 7.3)

(a) Estimate the integral F = [pic]dx using the sample mean MC method and 4000 trials. Your program should also compute the standard deviation from

( = [pic]

where

= [pic][pic] and = [pic][pic]

(b) The standard deviation computed above does not decrease with increasing n. Another estimate of the reliability of the MC-estimated value of F is from a standard deviation of mean values obtained from multiple simulations. Repeat the simulation of part (a) 19 more times and average the resultant values of F. From these 20 mean values compute (m, the standard deviation of the means

(m = [pic] (1)

(m can be estimated from a single run with the approximate relation

(m = (/201/2 (2)

Compare the value of (m computed from equations (1) and (2).

(c) Analytically determine the integral in part (a) and compare the exact result to the simulated result. Do you feel ( or (m is a better estimator of the standard deviation of a simulation run.

Solution

% Set 6, problem 4

%

n=4000;

f=0;fs=0;

for i=1:n

x=rand(1);fu=x*exp(-x);

f=f+fu;fs=fs+fu*fu;

end

area=f/n;areas=fs/n;

stde=sqrt(areas-area*area);stdee=stde/sqrt(20);

fprintf('Area = %8.4f Sigma = %8.4f\n',area,stde)

fprintf('Estimate standard deviation of the mean for 20 run, Sigma = Sigma/sqrt(20) = %8.6f\n',stdee)

f=zeros(20,1);fs=f;

for i=1:n

x=rand(20,1);fu=x.*exp(-x);

f=f+fu;fs=fs+fu.*fu;

end

area=f/n;areas=fs/n;

stde=sqrt(areas-area.*area);

for ir=1:10

fprintf('Area = %8.4f Sigma = %8.4f, Area = %8.4f, Sigma = %8.4f\n',area(ir),stde(ir),area(ir+10),stde(ir+10))

end

ave=sum(area)/20;

stdem=sqrt(sum((area-ave).^2)/19);

fprintf('Standard deviation of the mean for 20 run, Sigma = %8.6f\n',stdem)

>> s6p4

Area = 0.2649 Sigma = 0.1040

Estimate standard deviation of the mean for 20 run, Sigma = Sigma/sqrt(20) = 0.023245

Area = 0.2646 Sigma = 0.1041, Area = 0.2655, Sigma = 0.1040

Area = 0.2637 Sigma = 0.1055, Area = 0.2641, Sigma = 0.1049

Area = 0.2650 Sigma = 0.1048, Area = 0.2648, Sigma = 0.1044

Area = 0.2639 Sigma = 0.1045, Area = 0.2631, Sigma = 0.1052

Area = 0.2643 Sigma = 0.1053, Area = 0.2624, Sigma = 0.1047

Area = 0.2602 Sigma = 0.1069, Area = 0.2622, Sigma = 0.1066

Area = 0.2633 Sigma = 0.1053, Area = 0.2615, Sigma = 0.1060

Area = 0.2642 Sigma = 0.1057, Area = 0.2634, Sigma = 0.1053

Area = 0.2621 Sigma = 0.1046, Area = 0.2622, Sigma = 0.1058

Area = 0.2624 Sigma = 0.1068, Area = 0.2635, Sigma = 0.1053

Mean area = 0.2633

Standard deviation of the mean for 20 run, Sigma = 0.001319

(c) The analytical solution is

F = [pic]dx = [pic]= 1 ( 2e-1 = 0.2642

(m is a better estimator of the standard deviation of a simulation run

(m = 0.001319, ( = 0.1040

5. (Rowley, R.L., P. 7.4) Using importance sampling, compute the integral F = [pic]dx using 500, 1000, 5000, and 10,000 trials. Use P (x) = ax1/2 for the probability distribution function.

Solution

Normalize the probability distribution function

[pic]= 1 ( a[pic] = [pic]a[pic] = 1 ( a = [pic]

P (x) = [pic]x1/2

Q(x) = [pic]dx = [pic] = x3/2 = r ( x = r2/3

F = [pic]dx = [pic][pic] = [pic][pic][pic] = [pic][pic][pic]

% Set 6, problem 5

%

nv=[500 1000 5000 10000];

nt=length(nv);

for k=1:nt

f=0;n=nv(k);

for i=1:n

x=rand(1)^(2/3);fu=sqrt(x)*exp(-x);

f=f+fu;

end

area=2*f/(3*n);

fprintf('Area = %8.4f for %g trials\n',area,n)

end

>> s6p5

Area = 0.2633 for 500 trials

Area = 0.2639 for 1000 trials

Area = 0.2644 for 5000 trials

Area = 0.2641 for 10000 trials

Actual area = 0.2642

6. (Rowley, R.L., P. 7.5) Compute the integral F = [pic]dx

(a) using a uniform random number generator and (b) using importance sampling with

P (x) = [pic](1 ( x2)

for the probability distribution function. Use the same number of trials in the two cases and compare the variances (the square of the standard deviation).

Solution

The analytical solution is

F = [pic]dx = [pic]sin[pic] = [pic] = 0.63662

Q(x) = [pic]dx = [pic] = [pic][pic] = [pic]x ( [pic]= r

For a uniform random number r, x must be obtained from the cubic equation

x3 ( 3x + 2r = 0, then F = [pic]dx = [pic][pic]

% Set 6, problem 6

%

f=0;n=5000;

for i=1:n

x=rand(1);fu=cos(pi*x/2);

f=f+fu;

end

area=f/n;

fprintf('Uniform random number generator, Area = %8.4f\n',area)

f=0;

for i=1:n

r=rand(1); x=r;

for j=1:20

dx=(x^3-3*x+2*r)/(3*x*x-3);

x=x-dx;

if abs(dx)> s6p6

Uniform random number generator, Area = 0.6367

Importance sampling, Area = 0.6366

7. (Rowley, R.L., P. 7.6) Consider the object shown below which has been machined from a piece of material of uniform density. Find the center of mass of the object.

[pic]

Solution

We will choose random numbers for both x and y between ( 12 and 12. The density will be taken as zero if the coordinates place the point in one of the two holes in the disk or outside of the circle; otherwise the density is taken as one.

The mass of the disk per unit thickness is

M = ((122 ( 42 ( 22) = 389.56

The center of mass is given by

X = [pic][pic] = [pic][pic]

Y = [pic][pic] = [pic][pic]

We will generate random x and y and choose the density according to the algorithm:

If (x2 + y2) > 144 then ( = 0 because it is outside the circle

If ((x ( 4)2 + y2) < 16 then ( = 0 because it is in large hole

If (x2 + (y + 6) 2) < 4 then ( = 0 because it is in small hole

Otherwise ( = 1

% Set 6, problem 7

%

f=0;n=5000;

X=0;Y=0;

M=pi*(12^2-4^2-2^2);

for i=1:n

x=24*rand(1)-12;y=24*rand(1)-12;

den=1;

if (x*x+y*y)>144, den=0;end

if ((x-4)^2+y*y) s6p7

n = 5000 trials

X = -0.6710, Y = 0.1495

2 Rowley, R.L., Statistical Mechanics for Thermophysical Property Calculations, Prentice Hall, 1994

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download