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

### 17074190作业7

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.

(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

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

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

(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

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

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