EGR 511 - California State Polytechnic University, Pomona
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.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related searches
- california state university system
- california state university second bachelor s
- california state university tuition
- california state university jobs
- california state university system schools
- california state university system wiki
- california state university application log in
- california state university campuses list
- california state university log in
- california state university application deadline
- california state university tuition fee
- california state university fees