数值分析作业-三次样条插值

实验4.3三次样条插值函数（P126） 实验名称 实验时间 4.5三次样条插值函数的收敛性（P127） 姓名 班级 学号 成绩 实验4.3 三次样条差值函数

实验目的：

掌握三次样条插值函数的三弯矩方法。

实验函数:

f(x)12xet22dt

0.2 0.5793 0.3 0.6179 0.4 0.7554 x 0.0 0.1 F(x) 0.5000 0.5398 求f(0.13)和f(0.36)的近似值 实验内容：

（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；

（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线

比较插值结果。

实验4.5 三次样条差值函数的收敛性

实验目的：

多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。

实验内容：

按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。

实验要求：

（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情

况，分析所得结果并与拉格朗日插值多项式比较；

（2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考

xk yk  yk虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一

段数据如下： 0 1 2 3 4 5 6 7 8 9 10 0.0 0.8 0.79 1.53 2.19 2.71 3.03 3.27 2. 3.06 3.19 3.29 0.2 算法描述：

拉格朗日插值：

错误!未找到引用源。

其中错误!未找到引用源。是拉格朗日基函数，其表达式为：li(x)j0jinxxj(xixj)

牛顿插值：

Nn(x)f(x0)f[x0,x1](xx0)f[x0,x1,x2](xx0)(xx1)....

f[x0,x1,x2,...xn](xx0)(xx1)...(xxn1)f(xi)f(xj)f[xi,xj]xixjf[xj,xk]f[xi,xj]f[xi,xj,xk]xkxi其中.

.f[x0,x1...xn](f[x1,x2,...xn]f[x0,x1,...xn1])/(xnx0)三样条插值：

所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(aS(x)Mi1(xix)3(xxi1)3yyi1hi(MiMi1)Mi[i]x6hi6hihi6yhMhMy(i1ii1)(iii),x[xi1,xi]hi66hi

式中Mi=S(xi).

因此，只要确定了Mi的值，就确定了整个表达式，Mi的计算方法如下：

hiihihi1hi1令i

hihi1yi1yiyiyi16d()6f[xi1,xi,xi1]ihihi1hi1hi则Mi满足如下n-1个方程：

iMi12MiiMi1di,i1,2,...n1 常用的边界条件有如下几类：

m0,S(xn)ynmn （1） 给定区间两端点的斜率m0,mn,即S(x0)y0M0,S(xn)ynMn （2） 给定区间两端点的二阶导数M0，Mn,即S(x0)y0（3） 假设y=f(x)是以b-a为周期的周期函数，则要求三次样条插值函数S（x）

也为周期函数，对S（x）加上周期条件S(p)(x00)S(p)(xn0),p0,1,2

6y1y02MM(m0)01h1hi对于第一类边界条件有

M2M6(mnynyn1)n1nhnhn2M00M1d0对于第二类边界条件有

M2Mdnnnn1d060(f[x0,x1]m0)2(10)M0h16n(m0f[xn1,xn])2(1un)Mnhn其中

dn

那么解就可以为

2..0......M0d01.2..1.....M1d1.................M2d2 ......2..n1n1...............2nMn1dn1.Mdnn对于第三类边界条件，y0yn,M0Mn,S(x00)S(xn0)，由此推得

2M00M10Mn1d0，其中

0hnh16,0,d0(f[x0,x1]f[xn1,xn])，那么解就可以为： h1hnh1hnh1hn,,0.M0d02..0.......,Md.2.......1111..........M2d2....... ....n2..2..n2..Mdn2n..2n2n1..........dn1Mn1程序代码： 1拉格朗日插值函数

Lang.m

function f=lang(X,Y,xi) %X为已知数据的横坐标 %Y为已知数据的纵坐标 %xi插值点处的横坐标

%f求得的拉格朗日插值多项式的值 n=length(X); f=0; for i=1:n l=1; for j=1:i-1

l=l.*(xi-X(j))/(X(i)-X(j)); end; for j=i+1:n

l=l.*(xi-X(j))/(X(i)-X(j)); end;%拉格朗日基函数 f=f+l*Y(i); end

fprintf('%d\\n',f) return

2 牛顿插值函数

newton.m

function f=newton(X,Y,xi) %X为已知数据的横坐标 %Y为已知数据的纵坐标 %xi插值点处的横坐标

%f求得的拉格朗日插值多项式的值 n=length(X);

newt=[X',Y']; %计算差商表 for j=2:n for i=n:-1:1 if i>=j

Y(i)=(Y(i)-Y(i-1))/(X(i)-X(i-j+1)); else Y(i)=0; end end

newt=[newt,Y']; end

%计算牛顿插值 f=newt(1,2); for i=2:n z=1; for k=1:i-1 z=(xi-X(k))*z; end

f=f+newt(i-1,i)*z; end

fprintf('%d\\n',f) return

3三次样条插值第一类边界条件 Threch.m

function S=Threch1(X,Y,dy0,dyn,xi) % X为已知数据的横坐标 %Y为已知数据的纵坐标 %xi插值点处的横坐标 %S求得的三次样条插值函数的值 %dy0左端点处的一阶导数 % dyn右端点处的一阶导数 n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2);

for i=1:n%求函数的一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end

for i=2:n%求函数的二阶差商

f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end

d(1)=6*(f1(1)-dy0)/h(1);

d(n+1)=6*(dyn-f1(n-1))/h(n-1);%¸赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1

B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\\d; syms x; for i=1:n

Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4);

Sx(i)=vpa(Sx(i));%三样条插值函数表达式 end for i=1:n disp('S(x)=');

fprintf('%s (%d,%d)\\n',char(Sx(i)),X(i),X(i+1)); end for i=1:n

if xi>=X(i)&&xi<=X(i+1)

S=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3; end end

disp('xi S');

fprintf('%d,%d\\n',xi,S); return

4 三次样条插值第二类边界条件

Threch2.m

function [Sx]=Threch2(X,Y,d2y0,d2yn,xi)

X为已知数据的横坐标 %Y为已知数据的纵坐标 %xi插值点处的横坐标 %S求得的三次样条插值函数的值 %d2y0左端点处的二阶导数 % d2yn右端点处的二阶导数 n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2); for i=1:n%求一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end

for i=2:n%求二阶差商

f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end

d(1)=2*d2y0;

d(n+1)=2*d2yn;%赋初值 A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1

B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=0; A(n+1,n)=0; for i=1:n+1 A(i,i)=2; end for i=2:n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\\d; syms x; for i=1:n

Sx(i)=collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x-X(i))... +M(i)/2*(x-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x-X(i))^3); digits(4); Sx(i)=vpa(Sx(i));

end for i=1:n disp('S(x)=');

fprintf('%s (%d,%d)\\n',char(Sx(i)),X(i),X(i+1)); end for i=1:n

if xi>=X(i)&&xi<=X(i+1)

S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(xi-X(i))+M(i)/2*(xi-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(xi-X(i))^3; end end

disp('xi S');

fprintf('%d,%d\\n',xi,S); return

5插值节点处的插值结果 main3.m

clear clc

X=[0.0,0.1,0.2,0.3,0.4];

Y=[0.5000,0.5398,0.5793,0.6179,0.7554]; xi=0.13; %xi=0.36; disp('xi=0.13'); %disp('xi=0.36'); disp('拉格朗日插值结果'); lang(X,Y,xi); disp('牛顿插值结果'); newton(X,Y,xi);

disp('三次样条第一类边界条件插值结果');

Threch1(X,Y,0.40,0.36,xi);%0.4,0.36分别为两端点处的一阶导数 disp('三次样条第二类边界条件插值结果');

Threch2(X,Y,0,-0.136,xi);%0,-0.136分别为两端点处的二阶导数

6将多种插值函数即原函数图像画在同一张图上 main2.m

clear clc

X=[0.0,0.1,0.2,0.3,0.4];

Y=[0.5000,0.5398,0.5793,0.6179,0.7554]; a=linspace(0,0.4,21);

NUM=21; L=zeros(1,NUM); N=zeros(1,NUM); S=zeros(1,NUM); B=zeros(1,NUM); for i=1:NUM xi=a(i);

L(i)=lang(X,Y,xi);% 拉格朗日插值 N(i)=newton(X,Y,xi);% 牛顿插值

B(i)=normcdf(xi,0,1);%原函数

S(i)=Threch1(X,Y,0.4,0.36,xi);%三次样条函数第一类边界条件 end

plot(a,B,'--r'); hold on; plot(a,L,'b'); hold on; plot(a,N,'r'); hold on; plot(a,S,'r+'); hold on;

legend('原函数','拉格朗日插值','牛顿插值','三次样条插值',2); hold off

7增加插值节点观察误差变化 main4.m

clear; clc; N=5;

%4.5第一问

Ini=zeros(1,1001); a=linspace(-1,1,1001); Ini=1./(1+25*a.^2); for i=1:3 %节点数量变化次数 N=2*N;

t=linspace(-1,1,N+1);%插值节点 ft=1./(1+25*t.^2);%插值节点函数值 val=linspace(-1,1,101); for j=1:101

L(j)=lang(t,ft,val(j));

S(j)=Threch1(t,ft,0.074,-0.074,val(j));%三样条第一类边界条件插值 end

plot(a,Ini,'k')%原函数图象 hold on

plot(val,L,'r')%拉格朗日插值函数图像 hold on

plot(val,S,'b')%三次样条插值函数图像 str=sprintf('插值节点为%d时的插值效果',N); title(str);

legend('原函数','拉格朗日插值','三次样条插值');%显示图例 hold off figure end

8车门曲线 main5.m

clear clc

X=[0,1,2,3,4,5,6,7,8,9,10];

Y=[0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.,3.06,3.19,3.29]; dy0=0.8; dyn=0.2; n=length(X)-1; d=zeros(n+1,1); h=zeros(1,n-1); f1=zeros(1,n-1); f2=zeros(1,n-2);

for i=1:nh(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end

for i=2:nf2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end

d(1)=6*(f1(1)-dy0)/h(1);

d(n+1)=6*(dyn-f1(n-1))/h(n-1); A=zeros(n+1,n+1); B=zeros(1,n-1); C=zeros(1,n-1); for i=1:n-1

B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; for i=1:n+1 A(i,i)=2; end for i=2:n

A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\\d; x=zeros(1,n); S=zeros(1,n); for i=1:n x(i)=X(i)+0.5;

S(i)=Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i))*(x(i)-X(i))+M(i)/2*(x(i)-X(i))^2+(M(i+1)-M(i))/(6*h(i))*(x(i)-X(i))^3; end

plot(X,Y,'k'); hold on; plot(x,S,'o');

title('三次样条插值效果图');

legend('已知插值节点','三次样条插值'); hold off

实验结果：

4.3

1计算插值节点处的函数值

xi=0.13时

Xi=0.36时

2将多种插值函数即原函数图像画在同一张图上

4.5.1增加插值节点观察误差变化

从上面三张图可以看出增加插值节点并不能改善差之效果

4.5.2 车门曲线