﻿ 17074223-陈涛涛-作业七： 3.4and4.2Exercises - 计算数学达人 - 专，学者，数值代数，微分方程数值解

### 17074223-陈涛涛-作业七： 3.4and4.2Exercises

3.4 Exercises

1.Decide whether the equations form a cubic spline.

(b)

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

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

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

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

S1(x)’=12x+2

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

Because it satises all cubic spline properties for the data points (0, 1),so it can form a cubic spline.

(1, 2)

7.Solve equations (3.24) to nd the natural cubic spline through the three points (a) (0,0), (1,1), (2,4) (b) (1,1), (1,1), (2,4).

solution:

a)c1=0;c2=2/3;c3=0

b1=1/2;b2=2;

d1=1/2;d2=-1/2;

S1(x)=(1/2)x+(1/2)x^3

S2(x)=1+2(x-1)+(3/2)*(x-1)^2-(1/2)*(x-1)^3

b)c1=0;c2=2/3;c3=0

b1=-1;b2=2;

d1=1/4;d2=-1/2;

S1(x)=1-(x+1)+(1/4)*(x+1)^3

S2(x)=1+2(x-1)+(3/2)*(x-1)^2-(1/2)*(x-1)^3

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

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(n,n)=2;r(n)=vn;

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

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

%A(1,1:2)=[1 -1];  % parabol-term conditions, for n>=3

%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;       % solve for c coefficients

for i=1:n-1    % solve for b and d

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; % evaluate using nested multiplication

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)

4.2

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

a)t:0 1/4 1/2 3/4 y:1 3 2 0

b)t:0 1/4 1/2 3/4 y:1 3 2 1

c)t:0 1/2   1   3/2 y:3 1 3 2

solution:

a)A=[1 1 0;1 0 1;1 -1 0;1 0 -1];b=[1;3;2;0]

C=[3/2;-1/2;3/2]

y=3/2-1/2 cos 2πt + 3/2sin 2πt.

r=[0;0;0;0]

||r||2=0 ;RMSE=0;

b)A=[1 1 0;1 0 1;1 -1 0;1 0 -1];b=[1;3;2;1]

C=[7/4;-1/2;3/2]

y=7/4-1/2 cos 2πt + 3/2sin 2πt.

r=[-1/4;1/4;1/4;1/4]

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

c):A=[1 1 0;1 -1 0 ;1  1 0;1 -1 0];b=[3;1;3;2]

C=[9/4;3/4;3]

y=9/4-3/4 cos 2πt + 3sin 2πt.

r=[0;1/2;0;1/2]

||r||2=sqrt1/2 ;RMSE=sqrt(1/8);

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^c2ti .

a)t:-2 0 1 2 y:1 2 2 5

b)t:0 1 1 2 y:1 1 2 4

(a) y = 1.932e0.3615t, ||e||2 = 1.2825,

(b) y = 2t 1/4, ||e||2 = 0.9982

computer programs

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:

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 [c,RMSE]=lsquares(A,y)

a=A'*A;

b=A'*y;

c=gaussian(a,b);

r=y-A*c';

sum=0;

for i=1:length(r)

sum=sum+r(i)^2;

end

RMSE=sqrt(sum/length(r));

function x=gaussian(a,b)

n=length(b);

for j=1:n-1

if abs(a(j,j))<eps;error('zero pivot encountered');end

for i=j+1:n

mult=a(i,j)/a(j,j);

for k=j+1:n

a(i,k)=a(i,k)-mult*a(j,k);

end

b(i)=b(i)-mult*b(j);

end

end

for i=n:-1:1

for j=i+1:n

b(i)=b(i)-a(i,j)*x(j);

end

x(i)=b(i)/a(i,i);

end

A=[1 cos(pi/6)  sin(pi/6)  cos(pi/3);

1 cos(pi/3)  sin(pi/3)  cos(2*pi/3);

1 cos(pi/2)  sin(pi/2)  cos(pi);

1 cos(2*pi/3)  sin(2*pi/3)  cos(4*pi/3);

1 cos(5*pi/6)  sin(5*pi/6)  cos(5*pi/3);

1 cos(pi)      sin(pi)  cos(2*pi);

1 cos(7*pi/6)  sin(7*pi/6)  cos(7*pi/3);

1 cos(4*pi/3)  sin(4*pi/3)  cos(8*pi/3);

1 cos(3*pi/2)  sin(3*pi/2)  cos(3*pi);

1 cos(5*pi/3)  sin(5*pi/3)  cos(10*pi/3);

1 cos(11*pi/6)  sin(11*pi/6)  cos(11*pi/3);

1 cos(2*pi)  sin(2*pi)  cos(4*pi)];

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

>> [c,RMSE]=lsquares(A,y)

c =

5.5837    0.5921    0.4827   -0.0212

RMSE =

0.2284

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

>> A=[1 1960;1 1970;1 1990;1 2000];

>> y=[3039585530;3707475887;5281653820;6079603571];

>> [c,RMSE]=lsquares(A,log(y))

c =

-12.2623    0.0174

RMSE =

0.0147

>> format long

>> y=exp(1)^-12.2623*exp(1)^(0.0174*1980)

y =

4.333754095553840e+09

>> abs(y-4452584592)

ans =

1.188304964461603e+08