现在时间是:
当前位置:首 页 >> 数值分析(英)>> 教学区>> 文章列表

17074190作业7

作者:17074190   发布时间:2019-05-26 21:06:32   浏览次数:121

Exercise3.4

1.     问题:Decide whether the equations form a cubic spline.

(b) S(x)=  2x^3 + x^2 + 4x + 5, on [0,1]

(x 1)^3 + 7(x 1)^2 + 12(x 1) + 12, on [1,2]

Sol: We will check all three properties.

Property 1. There are n=3 data points. We must check

S1(0)=12 and S1(1)=12

S2(1)=12 and S2(2)=32  

Property 2. The rst derivatives of the spline functions are

S 1(x)=6*x^2+2*x+4

S 2(x)=3*(x-1)^2+14*(x-1)+12

We must check S 1(1)=S 2(1). The rst is 6+2+4= 12 ,which check out.

Property 3. The second derivatives are

S′′ 1(x)=12*x+2

S′′ 2(x)=6*(x-1)+14

We must check S′′ 1(1)=S′′ 2(1),which are true.Therefore,(b) is a cubic spline.

问题:7. Solve equations (3.24) to find the natural cubic spline through the three points

(a) (0,0),(1,1),(2,4)    (b) (-1,1),(1,1),(2,4)

Sol:

(a) The x-coordinates are x1=0,x2=1,and x3=2;The y-coordinates are a1=y1=0,

a2=y2=1,and a3=y3=4,and the differences are δ1=δ2=1,1=1,and 2=3.

The tridiagonal matrix equation (3.24) is

[ 1  0  0][c1]  [0]

[ 1  4  1][c2] = [6] ,the solution is [c1,c2,c3]=[0,3/2,0].Now,we get

[ 0  0  1][c3]  [0]

d1=(c2-c1)/3*δ1=1/2

d2=(c3-c2)/3*δ2= -1/2

b1=1/δ1- δ1(2*c1+c2)/3=1/2

b2=2/δ2- δ2(2*c2+c3)/3=2.

Therefore,the cubic spline is

S1(x)=0+1/2*x+0*x^2+1/2*x^3  on[0,1]

S2(x)=1+2*(x-1)+3/2*(x-1)^2-1/2*(x-1)^3 on[1,2]

(b)The x-coordinates are x1=-1,x2=1,and x3=2;The y-coordinates are a1=y1=1,

a2=y2=1,and a3=y3=4,and the differences are δ1=2,δ2=1,1=0,and 2=3.

The tridiagonal matrix equation (3.24) is

[ 1  0  0][c1]  [0]

[ 1  6  1][c2] = [9] ,the solution is [c1,c2,c3]=[0,3/2,0].Now,we get

[ 0  0  1][c3]  [0]

d1=(c2-c1)/3*δ1=1/4

d2=(c3-c2)/3*δ2= -1/2

b1=1/δ1- δ1(2*c1+c2)/3= -1

b2=2/δ2- δ2(2*c2+c3)/3=2.

Therefore,the cubic spline is

S1(x)=1-(x+1)+0*(x+1)^2+1/4*(x+1)^3  on[-1,1]

S2(x)=1+2*(x-1)+3/2*(x-1)^2-1/2*(x-1)^3 on[1,2]

3.4 Computer Problems

问题:1.Find the equations and plot the natural cubic spline that interpolation the data points (a) (0,3),(1,5),(2,4),(3,1)  (b) (-1,3),(0,5),(3,1),(4,1),(5,1).

Code:

function coeff=splinecoeff(x,y)

n=length(x);v1=0;vn=0;

A=zeros(n,n);

r=zeros(n,1);

for i=1:n-1

   dx(i)=x(i+1)-x(i);dy(i)=y(i+1)-y(i);

end

for i=2:n-1

   A(i,i-1:i+1)=[dx(i-1) 2*(dx(i-1)+dx(i)) dx(i)];

   r(i)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));

end

A(1,1)=1;

A(n,n)=1;

coeff=zeros(n,3);

coeff(:,2)=Ar;

for i=1:n-1

   coeff(i,3)=(coeff(i+1,2)-coeff(i,2))/(3*dx(i));

   coeff(i,1)=dy(i)/dx(i)-dx(i)*(2*coeff(i,2)+coeff(i+1,2))/3;

end

coeff=coeff(1:n-1,1:3);

function [x1,y1]=splineplot(x,y,k)

n=length(x);

coeff=splinecoeff(x,y);

x1=[];y1=[];

for i=1:n-1

   xs=linspace(x(i),x(i+1),k+1);

   dx=xs-x(i);

   ys=coeff(i,3)*dx;

   ys=(ys+coeff(i,2)).*dx;

   ys=(ys+coeff(i,1)).*dx+y(i);

   x1=[x1;xs(1:k)'];y1=[y1;ys(1:k)'];

end

x1=[x1;x(end)];y1=[y1;y(end)];

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

(a)

>> x1=[0 1 2 3];

>> y1=[3 5 4 1];

>> coeff=splinecoeff(x1,y1)

coeff =

   2.6667         0   -0.6667

   0.6667   -2.0000    0.3333

  -2.3333   -1.0000    0.3333

>> [x1,y1]=splineplot(x,y,4)

x1 =

        0

   0.2500

   0.5000

   0.7500

   1.0000

   1.2500

   1.5000

   1.7500

   2.0000

   2.2500

   2.5000

   2.7500

   3.0000

y1 =

   3.0000

   3.6563

   4.2500

   4.7188

   5.0000

   5.0469

   4.8750

   4.5156

   4.0000

   3.3594

   2.6250

   1.8281

1.0000

(b)

>> x1=[-1 0 3 4 5];

>> y1=[3 5 1 1 1];

>> coeff=splinecoeff(x1,y1)

coeff =

   2.5629         0   -0.5629

   0.8742   -1.6887    0.3176

  -0.6824    1.1698   -0.4874

   0.1950   -0.2925    0.0975

>> [x1,y1]=splineplot(x,y,5)

x1 =

        0

   0.2000

   0.4000

   0.6000

   0.8000

   1.0000

   1.2000

   1.4000

   1.6000

   1.8000

   2.0000

   2.2000

   2.4000

   2.6000

   2.8000

   3.0000

y1 =

   3.0000

   3.5280

   4.0240

  4.4560

   4.7920

   5.0000

   5.0560

   4.9680

   4.7520

   4.4240

   4.0000

   3.4960

   2.9280

   2.3120

   1.6640

1.0000

4.2 Exercises

问题:1.Fit data to the periodic model y = F3(t) = c1 + c2*cos2πt + c3*sin 2πt. Find the 2-norm error and the RMSE.

(a)

t y 0 1 1/4 3 1/2 2 3/4 0

(b)

t y 0 1 1/4 3 1/2 2 3/4 1

(c)

t y 0 3 1/2 1 1 3 3/2 2

Sol:

aThe corresponding inconsistent matrix equation is Ax =b, where  A=[1 cos0 sin0; 1 cos (2/π) sin (2/π); 1 cosπ sinπ; 1 cos (3π/2) sin (3π/2)=[1 1 0; 1 0 1 ;1 1 0; 1 0 -1] and b =[1;3;2;0]. The normal equations A* A*c =A* b are [4 0 0;0 2 0;0 1 2]*[c1;c2;c3]=[6;1;3]

which are easily solved as c1 =3/2, c2 =1/2, and c3 =7/4. The best version of the model, in the sense of least squares, is y =3/21/2*cos2πt+7/4*sin2πt,r=b-A*x=[1;3;2;0]-[1;13/4;2;-1/4]=[0;-1/4;0;1/4], ||r||2=r1^2+r2^2+r3^2+r4^2=1/2*2,  with RMSE= (||r||2)/4=1/4*2.

(b)The corresponding inconsistent matrix equation is Ax =b, where A=[1 cos0 sin0; 1 cos (2/π) sin (2/π); 1 cosπ sinπ; 1 cos (3π/2) sin (3π/2)=[1 1 0; 1 0 1 ;1 1 0; 1 0 -1] and b =[1;3;2;1]. The normal equations A* A*c =A* b are [4 0 0;0 2 0;0 1 2]*[c1;c2;c3]=[7;-1;0]

which are easily solved as c1 =7/4, c2 =1/2, and c3 =7/6. The best version of the model, in the sense of least squares, is y =7/41/2*cos2πt+7/6*sin2πt,r=b-A*x=[1;3;2;1]-[5/4;70/24;9/4;7/12]=[-1/4;1/12;-1/4;5/12], ||r||2=r1^2+r2^2+r3^2+r4^2=11/6, with RMSE= (||r||2)/4=11/12.

(c)The corresponding inconsistent matrix equation is Ax =b, where  A=[1 cos0 sin0; 1 cos π sin π; 1 cos2*π sin2*π; 1 cos 3π sin 3π]=[1 1 0; 1 -1 0 ;1 1 0; 1 -1 0] and b =[3;1;3;2]. The normal equations A’* A*c =A’* b are [4 0 0;0 4 0;0 0 0]*[c1;c2;c3]=[9;3;0]

which are easily solved as c1 =9/4, c2 =3/4, and c3 . The best version of the model, in the sense of least squares, is y =9/4+3/4*cos2πt+c3*sin2πt,r=b-A*x=[3;1;3;2]-[3;3/2;3;3/2]=[0;-1/2;0;1/2], ||r||2=√r1^2+r2^2+r3^2+r4^2=1/√2,  with RMSE= (||r||2)/√4=1/2*√2.

问题:3. Fit data to the exponential model by using linerization.Find the 2-norm of the differences between the data points yi and the best model c1*e^(c2^ti).

(a)

t y 2 1 0 2 1 2 2 5

(b)

t y 0 1 1 1 1 2 2 4

Sol:

a) Linearizing the model gives lny =k + c2t.

 k+(-2)*c2=0

 k+0*c2=ln2

 K+1*c2=ln2

 K+2*c2=ln5

A^T=[1 1 1 1;-2 0 1 2].

Because of A^T*A*x=A^T*b,we can get

[4 1;1 9][k;c2]=[ln20;ln50],

c2=0.36,k=0.66.

r=A*x-b=[0.06;0.66;1.02;1.38]-[0;0.69;0.69;1.6]=[-0.06;0.03;-0.33;0.22]

||r||2=(r1^2+r2^2+r^3+r^4)^1/2=0.4022

b) Linearizing the model gives lny =k + c2t.

 k+0*c2=0

 k+1*c2=0

 k+1*c2=ln2

 k+2*c2=ln4

A^T=[1 1 1 1;0 1 1 2].

Because of A^T*A*x=A^T*b,we can get

[4 4;4 6][k;c2]=[ln8;ln32]

C2=0.69,k=-0.17.

r=A*x-b=[0;0;ln2;ln4]-[0.17;0.52;0.52;1.21]=[0.17;-0.52;0,17;0.18]

||r||2=(r1^2+r2^2+r^3+r^4)^1/2=0.393

4.2 Computer Problems    

问题:1.Fit the monthly data for Japan 2003 oil consumption, shown in the following table, with the periodic model (4.9), and calculate the RMSE:

month   oil use (106 bbl/day)

Jan         6.224

Feb         6.665

Mar         6.241

Apr         5.302

May        5.073

Jun         5.127

Jul         4.994

Aug        5.012

Sep        5.108

Oct        5.377

Nov        5.510

Dec        6.372

Code

function LS(b,t)

n=length(t);

A=zeros(n,4);

for i=1:n

   A(i,1)=1;

   A(i,2)=cos(2*pi*t(i));

   A(i,3)=sin(2*pi*t(i));

   A(i,4)=cos(4*pi*t(i));

end

C=(A'*A)(A'*b)

>> t=[1;2;3;4;5;6;7;8;9;10;11;12];

>> b=[6.224;6.665;6.241;5.302;5.073;5.127;4.994;5.012;5.108;5.377;5.510;6.372];

>> LS(b,t)

C =

  1.0e+14 *

      NaN

      NaN

        0

  -7.9096

问题:3.Consider the world population data of Computer Problem 3.1.1. Find the best exponential fit of the data points by using linearization. Estimate the 1980 population, and find the estimation error.

Code:

function Model(t,b)

n=length(t);

A=zeros(n,2);

for i=1:n

   A(i,1)=1;

   A(i,2)=log(t(i));

end

C=(A'*A)(A'*b)

>> t=[1960;1970;1990;2000];

>> b=[log(3039585530);log(3707475887);log(5281653820);log(6079603571)];

>> Model(t,b)

C =

-239.3955

  34.4616

>> exp(-239.3955)

ans =

 1.0761e-104

>> y=(1.0761e-104)*1980^34.4616

y =

  4.3673e+09








上一篇:没有了    下一篇:没有了

Copyright ©2020    计算数学达人 All Right Reserved.

技术支持:自助建站 | 领地网站建设 |短信接口 |燕窝 版权所有 © 2005-2020 lingw.net.粤ICP备16125321号 -5