Download Assignment 8 - A Population Model for Numerical Computations | MATH 451 and more Assignments Mathematics in PDF only on Docsity!
CSE/MATH 451.002: Numerical Computations, spring
Answers to home work 8
1 Linear non-polynomial least square method
a). We should minimize
φ(α, β) =
∑^ m
k=
(α sin xk + β cos xk − yk)^2.
The normal equations become
∂φ ∂α
∑^ m
k=
2(α sin xk + β cos xk − yk) sin xk = 0
∂φ ∂β
∑^ m
k=
2(α sin xk + β cos xk − yk ) cos xk = 0
Rearrange it in a nice way, we get
α
sin^2 xk + β
sin xk cos xk =
yk sin xk α
sin xk cos xk + β
cos^2 xk =
yk cos xk
here all the summations go over k from 0 to m.
b). The normal equations give the linear system [ ∑ sin^2 xk
∑ sin^ xk^ cos^ xk sin xk cos xk
cos^2 xk
]
[
α β
]
∑ yk^ sin^ xk yk cos xk.
After we put in the values, we get [
- 88806 − 0. 332655 − 0. 332655 1. 11194
]
[
α β
]
which gives the solution: α = 0. 3335 , β = 3. 000.
2 A population model
a). In order to get the expression in the form y = a 0 + a 1 x we would rewrite the equation into 1 B
= KF T 0 − KF T
When we compute the values for (^) B^1 we see that they become very small. In order to make the computation easier we would rather use 10
12 B as the^ y−^ variable. Therefore we get the equation: 1012 B
= 10^12 KF T 0 − 1012 KF T
We still have an equation in the form y = a 0 + a 1 x. To make the computation even simpler we chose to centerize the Ti values, by introducing a variable change T˜ = T − T¯ = T − 1830. Now we have the following equation:
1012 B
= 10^12 KF (T 0 − 1830) − 1012 KF T˜
which is an equation in the form:
y = a 0 + a 1 x
where
y =
B
x = T − 1830 (2) a 0 = 1012 KF (T 0 − 1830) (3) a 1 = − 1012 KF (4)
This gives us the new table x y − 180 1835 − 130 1605 − 80 1374 − 30 1104 20 854 70 622 90 545 110 436 130 333 Solving this with the standard least squre method we get
a 0 =
yi 9
a 1 =
yixi ∑ x^2 i
The equation (3) and (4) give us now
KF = 4. 86 · 10 −^12 T 0 = 2029
0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
2
regression curve
o = observed values
Figure 1: Regression curve ploted together with the observation values.
Now you need to compute the quadratic sum
[yi − y(xi)]^2. The following command:
y_app=exp(0.690 -1.033x-1.450log(1+x))
would give us yi(xi). We can therefore set up the table:
xi yi yi 0 1. 996 1. 994
. 2 1. 244 1. 245 . 4. 810. 810 . 6. 541. 543 . 8. 375. 372
- 260
Use the command:
diff=(y-y_app).^2; sum(diff)
gives us the finally that
[yi − y(xi )]^2 = 1. 79 · 10 −^5.
Next is my suggestion to how you could solve the normal equations in Matlab:
% the values in the table x=[0.0:0.2:1.0]’; y=[1.996 1.244 0.810 0.541 0.375 0.259]’;
% Make a vector of ln(y_i) and ln(x_i+1): lny=log(y);
ln_nyx=log(x+1);
% Compute the elements in the system of normal equations: Acoeff=y A(1,1)=length(x); A(1,2)=sum(x); A(1,3)=sum(ln_nyx); A(2,1)=A(1,2); A(2,2)=x’x; A(2,3)=ln_nyx’x; A(3,1)=A(1,3); A(3,2)=A(2,3); A(3,3)=ln_nyx’ln_nyx;
rhs=[sum(lny) lny’x lny’ln_nyx]’; % right hand side
coeff=A\rhs;
%Get the original coefficients a0=exp(koeff(1)) a1=koeff(2) a2=-koeff(3)
% Plot the regression curve together with the observed values xer=[0.0:0.01:1.0]; func=a0exp(xer .a1)./(1+xer).^a2;
plot(xer,func,’y’,x,y,’ro’) gtext(’regression curve’) %when you run the file you need to click % with the mouse where the text should stay. gtext(’o = the observed values’)
4 Matlab least square method
The Matlab code looks like:
% Find the least square approximation to the following data set: x = [0:0.1:1]; y = [0.7829,0.8052,0.5753,0.5201,0.3783,0.2923,0.1695,0.0842,0.0415,0.009,0];
a1 = polyfit(x,y,1), t = [0:0.01:1]; v1 = polyval(a1,t); plot(x,y,’ro’,t,v1,’b’); title(’Regression with first order polynomial’), pause
a2 = polyfit(x,y,2), v2 = polyval(a2,t); plot(x,y,’ro’,t,v2,’b’); title(’Regression with second order polynomial’), pause
a4 = polyfit(x,y,4), v4 = polyval(a4,t); plot(x,y,’ro’,t,v4,’b’); title(’Regression with fourth order polynomial’), pause
a8 = polyfit(x,y,8), v8 = polyval(a8,t); plot(x,y,’ro’,t,v8,’b’); title(’Regression with eighth order polynomial’), pause
plot(x,y,’ro’,t,v1,’y-’,t,v2,’b:’,t,v4,’g.’,t,v8,’m-.’) title(’order 1: y-, order 2: b:, order 4: g., order 8: m-.’) print -depsc l8_fit.eps
And we get the following figure:
