﻿ 17074198 作业七 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074198 作业七

3.4Exercise

1. Problem: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]

Solution:

Property1:S1(0)=5,S1(1)=12

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

Property2:S1(x)=6x^2+2x+4,S2(x)=3(x-1)^2+14(x-1)+12

S1(1)=12,S2(1)=12

S1(1)=S2(1)

Property3:S1’’(x)=12x+2,S2’’(x)=6(x-1)+14

S1’’(1)=14,S2’’(1)=14

S1(1)=S2(1)

So it satisfied Property 1,2,3.

The equations is a cubic spline.

7. Problem: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).

Solution:(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;1,4,1;0,0,1][c1 c2 c3]=[0 6 0].

The solution is [c1, c2, c3] = [0,3/2,0]. Now, (3.22) and (3.23) yield

d1=(c2-c1)/3δ1=1/2,d2=(c3-c2)/3δ2=-1/2,

b1=△1/δ1-δ1(2c1+c2)/3=1/2,b2=△2/δ2-δ2(2c2+c3)/3=2.

Therefore, the cubic spline is S1(x)=x/2+x^3/2 on [0,1]

S2(x)=1+2(x-1)+3(x-1)^2/2-(x-1)^3/2 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;2,6,1;0,0,1][c1 c2 c3]=[0 9 0].

The solution is [c1, c2, c3] = [0,3/2,0]. Now, (3.22) and (3.23) yield

d1=(c2-c1)/3δ1=1/4,d2=(c3-c2)/3δ2=-1/6,

b1=△1/δ1-δ1(2c1+c2)/3=-1/4,b2=△2/δ2-δ2(2c2+c3)/3=2.

Therefore, the cubic spline is  S1(x)=1-(x+1)/4+(x+1)^3/4 on [-1,1]

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

3.4Computer Problems

1.Problem:Find the equations and plot the natural cubic spline that interpolates 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)

Result:

(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.2Exercises

1. Problem: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

Solution:(a)

Substituting the data into the model results in the following overdetermined system

of linear equations:

c1 + c2 cos2π(0)+ c3 sin 2π(0)=1

c1 + c2 cos2π(1/4) + c3 sin 2π(1/4)=3

c1 + c2 cos2π(1/2)+ c3 sin 2π(1/2)=2

c1 + c2 cos2πt(3/4)+ c3 sin 2π(3/4)=0

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]

b=[1 3 2 0].

The normal equations AT *Ac = AT* b are  [2,1,0,1;1,2,1,0;0,1,2,1; 1,0,1,2][c1 c2 c3 c4]=[6;-1;3].

c1=1.5,c2=-0.5,c3=1.5.

r=b-A*c=[1;3;2;0]-[1,1,0;1,0,1;1,-1,0;1,0,-1][1.5;-0.5;1.5]=[0;0;0;0]

||r||2=sqrt(r1^2+r2^2+r3^2+r4^2)=0

RMSE=(||r||2)/sqrt(m)=0.

(b)

Substituting the data into the model results in the following overdetermined system of linear equations:

c1 + c2*cos2π(0) + c3*sin 2π(0)=1

c1 + c2*cos2π(1/4) + c3*sin 2π(1/4)=3

c1 + c2*cos2π(1/2) + c3*sin 2π(1/2)=2

c1 + c2*cos2π(3/4) + c3*sin 2π(3/4)=1

The corresponding inconsistent matrix equation is Ax=b,where

[1   cos0     sin0  ]   [1  1  0]      [1]

A=[1  cos1/ sin1/] = [1  0  1] and b=[3].

[1   cosπ    sinπ ]  [1  -1  0]      [2]

[1  cos3/ sin3/]  [1  0  -1]      [1]

The normal equations AAc=Ab are

[4  0  0][c1]  [7]

[0  2  0][c2]= [-1],the solutions are c1 = 1.75, c2 = −0.5,and c3 = 1.

[0  0  2][c3]  [2]

[1]  [1  1  0][1.75] [-1/4]

r=b-Ac=[3] - [1  0  1][-0.5]= [1/4],

[2]  [1  -1  0][1]   [-1/4]

[1]  [1  0  -1]     [1/4]

The 2-norm error is SE=r1^2+r2^2+r3^2+r4^2=1/4,and the RMSE=sqrt(SE/m)=1/4.

(c)

Substituting the data into the model results in the following overdetermined system of linear equations:

c1 + c2*cos2π(0) + c3*sin 2π(0)=3

c1 + c2*cos2π(1/2) + c3*sin 2π(1/2)=1

c1 + c2*cos2π(1) + c3*sin 2π(1)=3

c1 + c2*cos2π(3/2) + c3*sin 2π(3/2)=2

The corresponding inconsistent matrix equation is Ax=b,where

[1   cos0     sin0  ]  [1  1  0]       [3]

A=[1   cosπ    sinπ ] = [1  -1  0] and b=[1].

[1   cos2π   sin2π]  [1  1  0]       [3]

[1   cos3π   sin3π]  [1  -1  0]      [2]

The normal equations AAc=Ab are

[4  0  0][c1]  [9]

[0  4  0][c2]= [3],the solutions are c1 =2.25, c2 = 0.75,and c3 = NaN.

[0  0  0][c3]  [0]

r=b-Ac 无解

The 2-norm error is none,and the RMSE is none.

2. Problem:Fit data to the exponential model by using linearization. Find the 2-norm of the difference between the data points yi and the best model c1*e^(c2*ti).

(a) t  y            (b)  t  y

−2 1                0  1

0  2                1  1

1  2                1  2

2  5                2  4

Solution:

(a) ln(c1*e^(c2*ti))=ln(c1)+c2*ti,let k=ln(c1),then ln(yi)=k+c2*ti.

Substituting the data into the linearization model yields

ln(1)=k+c2*(-2)

ln(2)=k+c2*(0)

ln(2)=k+c2*(1)

ln(5)=k+c2*2.

A=[1,-2;1,0;1,1;1,2],b=[ln1;ln2;ln2;ln5]

The normal equations AT *Ac = AT* b are [4,1;1,9][k;c2]=[2.9957;3.9120]

k=0.6586,c2=0.3615.c1=e^k=1.9321.

r=b-Ax=[0.0644;0.0346;-0.3269;0.2279]

||r||2=sqrt(r1^2+r2^2+r3^2+r4^2)≈0.4052

(b)

Substituting the data into the linearization model yields

k+c2^(0)=ln(1)

k+c2^(1)=ln(1)

k+c2^(1)=ln(2)

k+c2^(2)=ln(4)

The matrix equation is Ax=b,where x=(k,c2),

[1  0]      [ln(1)]

A=[1  1],and b=[ln(1)].The normal equations A’Ax=A’b are [4  4] [k] =[2.0794].

[1  1]      [ln(2)]                              [4  6][c2] [3.4657]

[1  2]      [ln(4)]

Which has solution k≈-0.1733 ,and c2≈0.6931,leading to c1=exp(k)≈0.8409.

r=b-Ax=[0.1733;-0.5199;0.1733;0.1733],

||r||2=sqrt(r1^2+r2^2+r3^2+r4^2)≈0.6003.

4.2Computer Problems

1.Problem: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)

Result:>> 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];

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

>> LS(b,t)

C =

1.0e+14 *

NaN

NaN

0

-7.9096

3.Problem: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)

Result:

>> 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