数值分析报告作业三次样条插值Word文档格式.docx
《数值分析报告作业三次样条插值Word文档格式.docx》由会员分享，可在线阅读，更多相关《数值分析报告作业三次样条插值Word文档格式.docx（20页珍藏版）》请在冰点文库上搜索。
理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：
（1）随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；
（2）三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考虑如下例子：
某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：
1
2
3
4
5
6
7
8
9
10
0.79
1.53
2.19
2.71
3.03
3.27
2.89
3.06
3.19
3.29
0.8
算法描述：
拉格朗日插值：
其中是拉格朗日基函数，其表达式为：
牛顿插值：
其中
三样条插值：
所谓三次样条插值多项式Sn（x）是一种分段函数，它在节点Xi（a<
X0<
X1……<
Xn<
b）分成的每个小区间[xi-1,xi]上是三次多项式，其在此区间上的表达式如下：
式中Mi=
.
因此，只要确定了Mi的值，就确定了整个表达式，Mi的计算方法如下：
令
则Mi满足如下n-1个方程：
常用的边界条件有如下几类：
（1）给定区间两端点的斜率m0,mn,即
（2）给定区间两端点的二阶导数M0，Mn,即
（3）假设y=f（x）是以b-a为周期的周期函数，则要求三次样条插值函数S（x）也为周期函数，对S（x）加上周期条件
对于第一类边界条件有
对于第二类边界条件有
那么解就可以为
对于第三类边界条件，
，由此推得
，其中
，那么解就可以为：
程序代码：
1拉格朗日插值函数
Lang.m
functionf=lang（X,Y,xi）
%X为已知数据的横坐标
%Y为已知数据的纵坐标
%xi插值点处的横坐标
%f求得的拉格朗日插值多项式的值
n=length（X）;
f=0;
fori=1:
n
l=1;
forj=1:
i-1
l=l.*（xi-X（j））/（X（i）-X（j））;
end;
forj=i+1:
%拉格朗日基函数
f=f+l*Y（i）;
end
fprintf（'
%d\n'
f）
return
2牛顿插值函数
newton.m
functionf=newton（X,Y,xi）
n=length（X）;
newt=[X'
Y'
];
%计算差商表
forj=2:
fori=n:
-1:
ifi>
=j
Y（i）=（Y（i）-Y（i-1））/（X（i）-X（i-j+1））;
elseY（i）=0;
newt=[newt,Y'
end
%计算牛顿插值
f=newt（1,2）;
fori=2:
z=1;
fork=1:
z=（xi-X（k））*z;
f=f+newt（i-1,i）*z;
return
3三次样条插值第一类边界条件
Threch.m
functionS=Threch1（X,Y,dy0,dyn,xi）
%X为已知数据的横坐标
%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）;
fori=1:
n%求函数的一阶差商
h（i）=X（i+1）-X（i）;
f1（i）=（Y（i+1）-Y（i））/h（i）;
end
fori=2:
n%求函数的二阶差商
f2（i）=（f1（i）-f1（i-1））/（X（i+1）-X（i-1））;
d（i）=6*f2（i）;
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）;
n-1
B（i）=h（i）/（h（i）+h（i+1））;
C（i）=1-B（i）;
A（1,2）=1;
A（n+1,n）=1;
n+1
A（i,i）=2;
n
A（i,i-1）=B（i-1）;
A（i,i+1）=C（i-1）;
M=A\d;
symsx;
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））;
%三样条插值函数表达式
disp（'
S（x）='
）;
fprintf（'
%s（%d,%d）\n'
char（Sx（i））,X（i）,X（i+1））;
ifxi>
=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;
xiS'
%d,%d\n'
xi,S）;
return
4三次样条插值第二类边界条件
Threch2.m
function[Sx]=Threch2（X,Y,d2y0,d2yn,xi）
X为已知数据的横坐标
%d2y0左端点处的二阶导数
%d2yn右端点处的二阶导数
n%求一阶差商
n%求二阶差商
d
（1）=2*d2y0;
d（n+1）=2*d2yn;
%赋初值
A（1,2）=0;
A（n+1,n）=0;
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;
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;
xi=0.13'
%disp（'
xi=0.36'
拉格朗日插值结果'
lang（X,Y,xi）;
牛顿插值结果'
newton（X,Y,xi）;
三次样条第一类边界条件插值结果'
Threch1（X,Y,0.40,0.36,xi）;
%0.4,0.36分别为两端点处的一阶导数
三次样条第二类边界条件插值结果'
Threch2（X,Y,0,-0.136,xi）;
%0,-0.136分别为两端点处的二阶导数
6将多种插值函数即原函数图像画在同一图上
main2.m
a=linspace（0,0.4,21）;
NUM=21;
L=zeros（1,NUM）;
N=zeros（1,NUM）;
S=zeros（1,NUM）;
B=zeros（1,NUM）;
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）;
%三次样条函数第一类边界条件
plot（a,B,'
--r'
holdon;
plot（a,L,'
b'
plot（a,N,'
r'
plot（a,S,'
r+'
legend（'
原函数'
'
拉格朗日插值'
牛顿插值'
三次样条插值'
2）;
holdoff
7增加插值节点观察误差变化
main4.m
clear;
clc;
N=5;
%4.5第一问
Ini=zeros（1,1001）;
a=linspace（-1,1,1001）;
Ini=1./（1+25*a.^2）;
3%节点数量变化次数
N=2*N;
t=linspace（-1,1,N+1）;
%插值节点
ft=1./（1+25*t.^2）;
%插值节点函数值
val=linspace（-1,1,101）;
101
L（j）=lang（t,ft,val（j））;
S（j）=Threch1（t,ft,0.,-0.,val（j））;
%三样条第一类边界条件插值
plot（a,Ini,'
k'
）%原函数图象
holdon
plot（val,L,'
）%拉格朗日插值函数图像
plot（val,S,'
）%三次样条插值函数图像
str=sprintf（'
插值节点为%d时的插值效果'
N）;
title（str）;
%显示图例
holdoff
figure
8车门曲线
main5.m
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;
nh（i）=X（i+1）-X（i）;
nf2（i）=（f1（i）-f1（i-1））/（X（i+1）-X（i-1））;
A=zeros（n+1,n+1）;
x=zeros（1,n）;
S=zeros（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;
plot（X,Y,'
holdon;
plot（x,S,'
o'
title（'
三次样条插值效果图'
已知插值节点'
holdoff
实验结果：
4.3
1计算插值节点处的函数值
xi=0.13时
Xi=0.36时
2将多种插值函数即原函数图像画在同一图上
4.5.1增加插值节点观察误差变化
从上面三图可以看出增加插值节点并不能改善差之效果
4.5.2车门曲线
升级会员 