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

### 17074201 作业七

P176

Exercise

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

Solution :

It’s not a cubic spline. We can know that S1’(x)=6x^2+2x+4,S2’(x)=3x^2+20x-11,

S1’(1)=12,S2’(1)=12,which satisfies the condition of cubic spline.

However,S1’’(x)=12x+2,S2’’(x)=6x+20,S1’’(1)=14,S2’’(1)=26,From that S1’’(1)doesn’t equal S2’’(1),we can conclude that it’s not a cubic spline.

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

Solution :
(a).

x1=0,x2=1,x3=2. a1=y1=0,a2=y2=1,a3=y3=4.

The tridiagonal matrix equation (3.24) is   &1=&2=1. 1=1，△2=3

The tridiagonal matrix equation 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,

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

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

b2=2/&2-&2/3(2c2+c3)=2

Therefore, the cubic spline is

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

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

(b).

x1=-1,x2=1,x3=2. a1=y1=1,a2=y2=1,a3=y3=4.

The tridiagonal matrix equation (3.24) is   &1=2,&2=1. 1=0，△2=3

The tridiagonal matrix equation 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,

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

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

b2=2/&2-&2/3(2c2+c3)=2

Therefore, the cubic spline is

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

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

Computer problem

Problem 1.

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

Solution :

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)); % right-hand side

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

(a).

>> x=[0 1 2 3];

>> y=[3 5 4 1];

>> coeff=splinecoeff(x,y)

coeff =

2.6667         0   -0.6667

0.6667   -2.0000    0.3333

-2.3333   -1.0000    0.3333

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

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

(b).

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

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

>> coeff=splinecoeff(x,y)

coeff =

2.5629         0   -0.5629

0.8742   -1.6887    0.3176

-0.6824    1.1698   -0.4874

0.1950   -0.2925    0.0975

S1(x)=3+2.5629(x+1)+0(x+1)^2-0.5629(x+1)^3    on[-1,0]

S2(x)=5+0.8742x-1.6887x^2+0.3176x^3         on[0,3]

S3(x)=1-0.6824(x-3)+1.1698(x-3)^2-0.4874(x-3)^3 on[3,4]

S4(x)=1+0.1950(x-4)-0.2925(x-4)^2+0.0975(x-4)^3 on[4,5]

P209

Exercise

Problem 1.

Fit data to the periodic model y =F3(t)=c1 + c2cos2πt + c3sin2π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:

1. a)c1+c2*1+c3*0=1

c1+c2*0+c3*1=3

c1+c2*(-1)+c3*0=2

c1+c2*0+c3*(-1)=0

[1 1 0;1 0 1;1 -1 0;1 0 -1]*[c1;c2;c3]=[1;3;2;0]

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

Because of A^T*A*x=A^T*b,

[4 0 0;0 2 0;0 0 2]*[c1;c2;c3]=[6;-1;3]

c1=3/2,c2=-1/2,c3=3/2

so y=3/2-1/2cos2πt+3/2sin2πt.

R=b-A*x=[1;3;2;0]-[1 1 0;1 0 1;1 -1 0;1 0 -1][3/2;-1/2;3/2]=[0;0;0;0]

||r||2=0,RMSE=||r||2/m^1/2=0

1. b)[1 1 0;1 0 1;1 -1 0;1 0 -1][c1;c2;c3]=[1;3;2;1]

A^T*A*x=A^T*b

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

[4 0 0;0 2 0;0 0 2]*[c1;c2;c3]=[7;-1;2]

c1=7/4,c2=-1/2,c3=1

so y=7/4-1/2cos2πt+sin2πt.

R=b-A*x=[1;3;2;1]-[1 1 0;1 0 1;1 -1 0;1 0 -1][7/4;-1/2;1]=[-1/4;1/4;-1/4;1/4]

||r||2=1/2,RMSE=1/2/2=1/4

1. c)       c1+c2*1+c3*0=3

c1+c2*(-1)+c3*0=1

c1+c2*1+c3*0=3

c1+c2*(-1)+c3*0=2

because of A^T*A*x=A^T*b,

[1 1 0;1 -1 0;1 10;1 -1 0][c1;c2;c3]=[3;1;3;2]

[4 0 0;0 4 0;0 0 0][c1;c2;c3]=[9;3;0]

c1=9/4,c2=3/4,c3 can be any constant in R.

so y=9/4+3/4cos2πt+c3sin2πt.

R=b-A*x=[3;1;3;2]-[1 1 0;1 -1 0;1 1 0;1 -1 0][9/4;3/4;c3]=[0;-1/2;0;1/2]

||r||2=√2/2,RMSE=√2/2/2=√2/4.

Problem 3.

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

(a)

t y −2 1 0 2 1 2 2 5

(b)

t y 0 1 1 1 1 2 2 4

Solution:

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

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

Computer problem

Problem 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

Solution :

function [x,r]=normequation(A,b)

m=pinv(A'*A);

n=A'*b;

x=m*n;

r=b-A*x;

Result:

>> x=(1+(0:11)/1)';

>> a=cos(2*pi*x);

>> b=sin(2*pi*x);

>> c=cos(4*pi*x);

>> d=ones(12,1);

>> A=[d a b c]

A =

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

1.0000    1.0000   -0.0000    1.0000

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

b =

6.2440

6.6650

6.2410

5.3020

5.0730

5.1270

4.9940

5.0120

5.1080

5.3770

5.5100

6.3720

>> [x,r]=normequation(A,b)

x =

1.8618

1.8618

-0.0000

1.8618

r =

0.6586

1.0796

0.6556

-0.2834

-0.5124

-0.4584

-0.5914

-0.5734

-0.4774

-0.2084

-0.0754

0.7866

Problem 3.

Consider the world population data of Computer Problem 3.1.1. Find the best exponential ﬁt of the data points by using linearization. Estimate the 1980 population, and ﬁnd the estimation error.