﻿ 17074199 3.4+4.2 第7次作业 - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074199 3.4+4.2 第7次作业

P176  3.4 Exercises

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]

PROPERTY1   S1(0)=5     S1(1)=12

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

PROPERTY2   The first derivatives of the spline functions are

S1’(x)=6x^2+2x+4

S2’(x)=3(x-1)^2+14(x-1)+12

We must check S1’(1)= S2’(1).

S1’(1)= 12      S2’(1)=12

S1’(1)= S2’(1)

It satisfied with property 2.

PROPERTY3   The second derivatives of the spline functions are

S1’’(x)=12x+2

S2’’(x)=6(x-1)+14

We must check S1’’(1)= S2’’(1).

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

S1’’(1)= S2’’(1)

It satisfied with property 3.

Therefore, the equations 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).

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

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+x^3/2              on [0,1]

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

The solution is [c1;c2;c3] = [0;3/2;0]

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)+1/4(x+1)^3     on [-1,1]

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

3.4 Computer Problems

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

(a) (0,3), (1,5), (2,4), (3,1)

%Program 3.5 Calculation of spline coefficients

%Calculates coefficients of cubic spline

%Input:x,y vectors of data points

% plus two optional extra data v1,vn

%Output:matrix of coeffcients b1,c1,d1;b2,c2,d2;¡­

function coeff=splinecoeff(x,y)

n=length(x);

v1=0;

vn=0;

A=zeros(n,n);         %matrix A is nxn

r=zeros(n,1);

for i=1:n-1           %define the deltas

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

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

end

for i=2:n-1           %load the A matrix

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

%Set endpoint conditions

%Use only one of following 5 pairs:

A(1,1)=1;             %natural spline conditions

A(n,n)=1;

%A(1,1)=2;r(1)=v1;

%A(n,n)=2;r(n)=vn;

%A(1,1:2)=[2*dx(1) dx(1)]; r(1)=3*(vn-dy(n-1)/dx(n-1)); %clamped

%A(n,n-1:n)=[1 -1];

%A(1,1:3)=[dx(2) -(dx(1)+dx(2)) dx(1)]; %not-a-knot,for n>=4

%A£¨n,n-2:n)=[dx(n-1) -(dx(n-2)+dx(n-1)) dx(n-2)];

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)

>> x=[0 1 2 3];

y=[3 5 4 1];

coeff=splinecoeff(x,y)

coeff =

2.666666666666667                   0  -0.666666666666667

0.666666666666667  -2.000000000000000   0.333333333333333

-2.333333333333334  -1.000000000000000   0.333333333333333

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

x1 =

0

0.250000000000000

0.500000000000000

0.750000000000000

1.000000000000000

1.250000000000000

1.500000000000000

1.750000000000000

2.000000000000000

2.250000000000000

2.500000000000000

2.750000000000000

3.000000000000000

y1 =

3.000000000000000

3.656250000000000

4.250000000000000

4.718750000000000

5.000000000000000

5.046875000000000

4.875000000000000

4.515625000000000

4.000000000000000

3.359375000000000

2.625000000000000

1.828125000000000

1.000000000000000

(b) >> x=[-1 0 3 4 5]

x =

-1     0     3     4     5

>> y=[3 5 1 1 1]

y =

3     5     1     1     1

>> coeff=splinecoeff(x,y)

coeff =

2.562893081761006                   0  -0.562893081761006

0.874213836477987  -1.688679245283019   0.317610062893082

-0.682389937106918   1.169811320754717  -0.487421383647799

0.194968553459119  -0.292452830188679   0.097484276729560

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

x1 =

-1.000000000000000

-0.800000000000000

-0.600000000000000

-0.400000000000000

-0.200000000000000

0

0.600000000000000

1.200000000000000

1.800000000000000

2.400000000000000

3.000000000000000

3.200000000000000

3.400000000000000

3.600000000000000

3.800000000000000

4.000000000000000

4.200000000000000

4.400000000000000

4.600000000000000

4.800000000000000

5.000000000000000

y1 =

3.000000000000000

3.508075471698113

3.989132075471698

4.416150943396226

4.762113207547170

5.000000000000000

4.985207547169812

4.166188679245283

2.954566037735849

1.761962264150943

1.000000000000000

0.906415094339623

0.883018867924528

0.906415094339623

0.953207547169811

1.000000000000000

1.028075471698113

1.037433962264151

1.032754716981132

1.018716981132076

1.000000000000000

4.2 Exercises

1.Fit data to the periodic model y = F 3 (t) = c 1 + c 2 cos2πt + c 3 sin2πt. Find the 2-norm error

and the RMSE.

(a)    t=0 1/4 1/2 3/4

y=1  3  2  0

The system of equations is now

c1+c2* cos2π(0)+c3* sin2π(0)

c1+c2* cos2π(1/4)+c3* sin2π(1/4)

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

c1+c2* cos2π(3/4)+c3* sin2π(3/4)

The corresponding inconsistent matrix equation is Ax=b,where:

A=[1 cos0 sin0]=[1 1 0]

=[1 cos (π/2) sin(π/2)]=[1 0 1]

=[1 cos (π) sin(π)]=[1 -1 0]

=[ cos (3π/2) sin(3π/2)]=[1 0 -1]

And b=[1;3;2;0]

A^T *A*c=A^T*b leading to the following normal equations:

A^T *A=[1 1 1 1;1 0 -1 0;0 1 0 -1][1 1 0;1 0 1;1 -1 0;1 0 -1]=[4 0 0;0 2 0;0 0 2]

A^T*b=[1 1 1 1;1 0 -1 0;0 1 0 -1][1;3;2;0]=[6;-1;3]

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

C1=1.5

C2=-0.5

C3=1.5

t      y     line    error

0      1       1       0

1/4    3       3       0

1/2    2       2       0

3/4    0       0       0

The best version of the model,in the sense of least squares, is y=1.5-0.5 cos2πt+1.5 sin2πt

,with RMSEsqrt(SE/m)=sqrt((r1^2+r2^2+r3^2+r4^2)/4)=0

(b)    t=0 1/4 1/2 3/4

y=1  3  2  1

The system of equations is now

c1+c2* cos2π(0)+c3* sin2π(0)

c1+c2* cos2π(1/4)+c3* sin2π(1/4)

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

c1+c2* cos2π(3/4)+c3* sin2π(3/4)

The corresponding inconsistent matrix equation is Ax=b,where:

A=[1 cos0 sin0]=[1 1 0]

=[1 cos (π/2) sin(π/2)]=[1 0 1]

=[1 cos (π) sin(π)]=[1 -1 0]

=[ cos (3π/2) sin(3π/2)]=[1 0 -1]

And b=[1;3;2;1]

A^T *A*c=A^T*b leading to the following normal equations:

A^T *A=[1 1 1 1;1 0 -1 0;0 1 0 -1][1 1 0;1 0 1;1 -1 0;1 0 -1]=[4 0 0;0 2 0;0 0 2]

A^T*b=[1 1 1 1;1 0 -1 0;0 1 0 -1][1;3;2;1]=[7;-1;2]

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

C1=1.75

C2=-0.5

C3=1

t      y     line    error

0      1     1.25   -0.25

1/4    3      2.75   0.25

1/2    2      2.25   -0.25

3/4    1      0.75   0.25

The best version of the model,in the sense of least squares, is y=1.75-0.5 cos2πt+sin2πt

,with RMSEsqrt(SE/m)=sqrt((r1^2+r2^2+r3^2+r4^2)/4)=0.25

(c)    t=0 1/4 1/2 3/4

y=3  1  3  2

The system of equations is now

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 A’Ac=A’b 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]

The best version of the model,in the sense of least squares, is y=2.25+0.75*cos2πt

r=b-Ac 无解

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

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 c1*e^(c2*ti) .

(a)    t -2 0 1 2

y 1 2 2 5

(b)    t 0 1 1 2

y 1 1 2 4

(a)    y=c1*e^(c2*t)

Iny=k+c2*t

Substituting the data into the linearized model yields

k+c1(-2)=In(1)

k+c1(0)=In(2)

k+c1(1)=In(2)

k+c1(2)=In(5)

and so forth.The matrix equation is Ax=b,where x=(k,c2),

where we have set k=In(c1).This leads to the matrix equation Ax=b,where

A=[1 -2;1 0;1 1;1 2] and b=[In(1);In(2);In(2);In(5)]

The normal equations allow determination of k and c2,and c1=e^k.

[1  -2]      [ln(1)]

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

[1  1]      [ln(2)]                               [1  9][c2] [3.9120]

[1  2]      [ln(5)]

Which has solution k≈0.6586 and c2≈0.3615,leading to c1=exp(k)≈1.9321.

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

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

(b)    y=c1*e^(c2*t)

Iny=k+c2*t

Substituting the data into the linearized model yields

k+c1(-2)=In(1)

k+c1(0)=In(2)

k+c1(1)=In(2)

k+c1(2)=In(5)

and so forth.The matrix equation is Ax=b,where x=(k,c2),

where we have set k=In(c1).This leads to the matrix equation Ax=b,where

A=[1 -2;1 0;1 1;1 2] and b=[In(1);In(2);In(2);In(5)]

The normal equations allow determination of k and c2,and c1=e^k.

[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.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(10^6 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:

y=c1+c2cos2πt +c3* sin2πt+ c4cos4πt

y=c1+c2cos2π(1) +c3* sin2π(1)+ c4cos4π(1)

y=c1+c2cos2π(2) +c3* sin2π(2)+ c4cos4π(2)

y=c1+c2cos2π(3) +c3* sin2π(3)+ c4cos4π(3)

……

function LST(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];

LST(b,t)

C =

1.0e+14 *

NaN

NaN

0

-7.9096

r=b-A*C

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.

year  population

1960  3039585530

1970  3707475887

1990  5281653820

2000  6079603571

Code:

We use the the best exponential fit of the data points y=c1*t^(c2)

In the case of the exponential model, the model is lineared by applying the natural logarithm:

Iny=Inc1+c2*Int=k+c2*Int

A=[1 1960;1 1970;1 1990;1 2000] and b=[In3039585530;In3707475887;In5281653820;In6079603571]

Solving the function Ac=b

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 =

1.0e+02 *

0.344615755264695

>> exp(-2.393955462904132)

ans =

0.091267961922513

>> t1=1980

t1 =

1980

>> y= 0.091267961922513*t1.^(0.344615755264695)

y =

1.248536137189043

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

>> C=[-2.393955462904132; 0.344615755264695]

C =

-2.393955462904132

0.344615755264695

>> r=b-A*C

r =

16.704261929492191

16.834440983198466

17.066373590342820

17.158586297470215

Solving the linear least squares problem yields k=-2.393955462904132

c2= 0.344615755264695 c1= 0.091267961922513

The equation are:

Iny=Inc1+c2*Int=k+c2*Int

Iny=-2.393955462904132+0.344615755264695*Int

y= 0.091267961922513*t^(0.344615755264695)

We estimate the 1980 population is 1.248536137189043

The residual of the least squares solution c as

r=b-A*C=

16.704261929492191

16.834440983198466

17.066373590342820

17.158586297470215

However, I find the best exponential fit of the data points describing the polpulation is y=c1*e^(c2*t)

So we should use the the best exponential fit of the data points y=c1*e^(c2*t)

In the case of the exponential model, the model is lineared by applying the natural logarithm:

Iny=Inc1+c2*t=k+c2*t

A=[1 1960;1 1970;1 1990;1 2000] and b=[In3039585530;In3707475887;In5281653820;In6079603571]

Solving the function Ac=b

function Model(t,b)

n=length(t);

A=zeros(n,2);

for i=1:n

A(i,1)=1;

A(i,2)=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 =

-12.262349027671746

0.017403246298083

The following solution is the same as the above context.