数值分析作业-三次样条插值教学提纲

数值分析作业-三次样
条插值
数值计算方法作业
实验4.3 三次样条差值函数
实验目的：
掌握三次样条插值函数的三弯矩方法。

实验函数:
dt e
x f x
t ⎰
∞
--
=2
221)(π
求f(0.13)和f(0.36)的近似值
实验内容：
（1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值；
（3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲
线比较插值结果。

实验4.5 三次样条差值函数的收敛性
实验目的：
多项式插值不一定是收敛的，即插值的节点多，效果不一定好。

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

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

实验要求：
（1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情
况，分析所得结果并与拉格朗日插值多项式比较；
（2） 三次样条插值函数的思想最早产生于工业部门。

作为工业应用的例子，
考虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：
k x 0 1 2 3 4 5 6 7 8 9 10 k y 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29
k
y ' 0.8
0.2
算法描述：
拉格朗日插值：
其中
是拉格朗日基函数，其表达式为：()
∏
≠=--=n
i j j j i j
i x x x x x l 0
)
()(
牛顿插值：
)
)...()(](,...,,[....
))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N
其中⎪⎪⎪⎪⎪
⎪⎩
⎪⎪
⎪⎪⎪⎪⎨⎧
--=--=
--=
-)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i j
i j i j i
三样条插值：
所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点
Xi(a<X0<X1……<Xn<b)分成的每个小区间[x i-1,x i ]上是三次多项式，其在此区间上的表达式如下：
],[),6()6(
]6
)
([6)(6)()(111113131i i i
i i i i i i i i i i i i i i i i i i x x x h y
M h M h h y x M M h h y y h x x Mi h x x M x S -------∈-+-+---+-+-=
式中Mi=)(i x S ''.
因此，只要确定了Mi 的值，就确定了整个表达式，Mi 的计算方法如下：
令⎪⎪⎪⎪⎩
⎪⎪⎪⎪⎨⎧=---+=
+=+=+--++++++],,[6)(6
111111111
i i i i i i i i i i i i i i i i i i i i
x x x f h y y h y y h h d h h h h h h λμ
则Mi 满足如下n-1个方程：
1,...2,1,211-==+++-n i d M M M i i i i i i λμ 常用的边界条件有如下几类：
（1） 给定区间两端点的斜率m 0,m n ,即n n n m y x S m y x S ='='='=')(,)(00
（2） 给定区间两端点的二阶导数M0，Mn,即
n n n M y x S M y x S =''=''=''='')(,)(00
0 （3） 假设y=f(x)是以b-a 为周期的周期函数，则要求三次样条插值函数S
（x ）也为周期函数，对S （x ）加上周期条件
2,1,0),0()0()(0)(=-=+p x S x S n p p
对于第一类边界条件有⎪⎪⎩
⎪⎪⎨⎧
--=+--=+--)(62)(6211001110n n n
n n
n i
h y y mn h M M m h y y h M M
对于第二类边界条件有⎩⎨⎧=+=+-n
n n n d M M d M M 2210
100μλ
其中
n
n n n n
n n M u x x f m h d M m x x f h d )1(2]),[(6)1(2)],[(610000101
0-+-=-+-=
-μ
λλ
那么解就可以为
⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d d M M M M M 1210121011...2...............2............................1..2.1......0..2μλμλμλ 对于第三类边界条件，)0()0(,,000-=+==n n n x S x S M M y y ，由此推得
0010012d M M M n =-++μλ，其中
]),1[],[(6
,,101010110n n n
n n n x x f x x f h h d h h h h h h --+=+=+=
μλ，那么解就可以为：
⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢
⎢⎢
⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-------1221012101221100...2.............2..............................2..,,.......,
..22
n n n n n n n d d d d d M M M M M n μλλμλμμλ 程序代码： 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.89,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 车门曲线。