数值分析作业三次样条插值.docx
《数值分析作业三次样条插值.docx》由会员分享，可在线阅读，更多相关《数值分析作业三次样条插值.docx（23页珍藏版）》请在冰豆网上搜索。
数值分析作业三次样条插值
数值计算方法作业
实验名称
实验4.3三次样条插值函数
（P126）
4.5三次样条插值函数的收敛性
（P127）
实验时间
姓名
班级
学号
成绩
实验4.3三次样条差值函数
实验目的：
掌握三次样条插值函数的三弯矩方法
实验函数：
求f（0.13）和f（0.36）的近似值
实验内容:
（1）编程实现求三次样条插值函数的算法，分别考虑不同的边界条件；
（2）计算各插值节点的弯矩值；
（3）在同一坐标系中绘制函数f（x），插值多项式，三次样条插值多项式的曲线
比较插值结果。
实验4.5三次样条差值函数的收敛性
实验目的：
多项式插值不一定是收敛的，即插值的节点多，效果不一定好。
对三次样条插值函数如何呢？
理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。
实验内容：
按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。
实验要求：
（1）随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情况，分析所得结果并与拉格朗日插值多项式比较；
（2）三次样条插值函数的思想最早产生于工业部门。
作为工业应用的例子，考
虑如下例子：
某汽车制造商根据三次样条插值函数设计车门曲线，其中一段数据如下：
Xk
0
1
2
3
4
5
6
7
8
9
10
yk
0.0
0.79
1.53
2.19
2.71
3.03
3.27
2.89
3.06
3.19
3.29
yk
0.8
0.2
算法描述:
拉格朗日插值：
错误！
未找到引用源。
n（x_X）其中错误！
未找到引用源。
是拉格朗日基函数，其表达式为：
h（x）」
j=0（xi-Xj）
牛顿插值：
Nn（x）=f（Xg）f[Xo,Xi]（X-xO）f[Xo,Xi,X2〕（X-xO）（x-Xi）•…
hi
hiMi4hiMiyi
-6）（6*,皿"]
式中Mi=S（Xi）.
f[Xg,Xi...Xn]=（f[Xi,X2,...Xn]-f[X。
，为，..人」）/（X.-Xg）
三样条插值:
所谓三次样条插值多项式Sn（x）是一种分段函数，它在节点
Xi（a.Mi（x—Xy）3.[yi-yi4hi（Mi-My）
6hihi6
则Mi满足如下n-1个方程:
7Mi」■2Mi…冷Mii=di,i=1,2,...n—'1
常用的边界条件有如下几类：
（1）给定区间两端点的斜率mo,mn,即s（x0）=y0=m0,S（xn）=yn=mn
（2）给定区间两端点的二阶导数MO,Mn,即S（XcHy。
=M。
S”（Xn）=y^Mn
（3）假设y=f（X）是以b-a为周期的周期函数，则要求三次样条插值函数S（x）
也为周期函数，对S（x）加上周期条件S（p）（x0•0）=S（p）（Xn-0）,p=01,2
6y1—y0
2M0+M1=—-m0）
对于第一类边界条件有
h1hi
Mn」2Mn6（mn归）
hnhn
那么解就可以为
'2.7^0.....
田.2.21
Ym°
M1
M2
d0
d1
d2
n4..2..n』
Jn..2
Hn
丿LM
对于第三类边界条件,
y。
=yn,M°二Mn,S（X。
0）=S（Xn-0），由此推得
2M0…0M1」°MnT二d°，其中
hl,J0hn,do6（fix。
”]—f[Xn-1,Xn]），那么解就可以为:
hi-hnhihnhi-hn
■2•払…
•••,,,2
Mo1
'do1
已22
1……
M1
d1
M2
d2
LL…•厂n_2
•2・・'-n_2
=
入n」……
•••%.21
M
n_2_
dn/
）
'Mn—
1
dn/
程序代码：
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:
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
functionf=newton（X,Y,xi）
%X为已知数据的横坐标
%Y为已知数据的纵坐标
%xi插值点处的横坐标
%f求得的拉格朗日插值多项式的值
n=length（X）;
newt=[X',Y'];
%计算差商表
forj=2:
n
fori=n:
-1:
1
ifi>=j
丫（i）=（丫（i）-Y（i-1））/（X（i）-X（i-j+1））;
elseY（i）=0;
end
end
newt=[newt,Y'];
end
%计算牛顿插值
f=newt（1,2）;
fori=2:
n
z=1;
fork=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
functionS=Threch1（X,Y,dyO,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）;
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）;
end
d
（1）=6*（f1⑴-dyO）/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）;
fori=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;
fori=1:
n+1
A（i,i）=2;
end
fori=2:
n
A（i,i-1）=B（i-1）;
A（i,i+1）=C（i-1）;
end
M=A\d;
symsx;
fori=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））A2+（M（i+1）-M（i））/（6*h（i））*（x-X（i））A3）;
digits⑷；
Sx（i）=vpa（Sx（i））;%三样条插值函数表达式
end
fori=1:
n
disp（'S（x）='）；
fprintf（'%s（%d,%d）\n',char（Sx（i））,X（i）,X（i+1））;
end
fori=1:
n
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））A2+（M（i+1）-M（i））/（6*h（i））*（xi-X（i））A3;
end
end
disp（'xiS'）;
fprintf（'%d,%d\n',xi,S）;
return
4三次样条插值第二类边界条件
Threch2.m
function[Sx]=Threch2（X,Y,d2yO,d2yn,xi）
X为已知数据的横坐标
%Y为已知数据的纵坐标
%xi插值点处的横坐标
%殊得的三次样条插值函数的值
%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）;
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）;
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）;
fori=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;
fori=1:
n+1
A（i,i）=2;
end
fori=2:
n
A（i,i-1）=B（i-1）;
A（i,i+1）=C（i-1）;
end
M=A\d;
symsx;
fori=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））A2+（M（i+1）-M（i））/（6*h（i））*（x-X（i））A3）;
digits⑷；
Sx（i）=vpa（Sx（i））;
end
fori=1:
n
disp（'S（x）='）;
fprintf（'%s（%d,%d）\n',char（Sx（i））,X（i）,X（i+1））;
end
fori=1:
n
ifxi>=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）F2+（M（i+1）-M（i））/（6*h（i））*（xi-X（i））A3;
end
end
disp（'xiS'）;
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（'拉格朗日插值结果’）；
Iang（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）;
fori=1:
NUM
xi=a（i）;
%拉格朗日插值
%牛顿插值
%原函数
%三次样条函数第一类边界条件
L（i）=lang（X，丫,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'
）；
holdon;plot（a,L,
'b'
）；
holdon;
plot（a,N,
holdon;
'r'
）；
plot（a,S,
holdon;
legend（'原函数’，’拉格朗日插值’，‘牛顿插值’,'三次样条插值’,2）;
holdoff
'r+'
）;
7增加插值节点观察误差变化
main4.m
clear;
clc;
N=5;
%4.5第一问
Ini=zeros（1,1001）;
a=linspace（-1,1,1001）;
Ini=1./（1+25*a.A2）;
%节点数量变化次数
fori=1:
3
N=2*N;
t=linspace（-1,1,N+1）;
ft=1./（1+25*t.A2）;val=linspace（-1,1,101）;
forj=1:
101
L（j）=lang（t,ft,val（j））;
%插值节点
%插值节点函数值
S（j）=Threch1（t,ft,0.074,-0.074,val（j））;
end
%三样条第一类边界条件插值
plot（a,Ini,holdon
'k'）%原函数图象
plot（val,L,'r'）%拉格朗日插值函数图像
holdon
plot（val,S,'b'）%三次样条插值函数图像
str=sprintf（'插值节点为％d时的插值效果',N）;
title（str）;
legend（'原函数’，’拉格朗日插值’，‘三次样条插值’）；％显示图例
holdoff
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）;
fori=1:
nh（i）=X（i+1）-X（i）;
f1（i）=（Y（i+1）-Y（i））/h（i）;
end
fori=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）;
fori=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;
fori=1:
n+1
A（i,i）=2;
end
fori=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）;
fori=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））A2+（M（i+1）-M（i））/（6*h（i））*（x（i）-X（i））A3;
end
plot（X,Y,'k'）;holdon;
plot（x,S,'o'）；
title（'三次样条插值效果图’）；
legend（'已知插值节点’，‘三次样条插值’）;
holdoff
实验结果：
4.3
1计算插值节点处的函数值
xi=0.13时
k±=0-13
才立倍钢曰插危镭果
5.S36&03e-301
半硕晒恺结果
5,665517e-aOl
三样条弟——榮边界寿^牛価值蜡果
.5000+.400a*K-.3S98*x^2+3.698*k*3（O,1,OOOOOOe-OO1>
S>=
-514T-1O,93*k^3-i-4,0IE冰x"2一”<1”OOOOOOe-OOi2”ODOOOOe-OO1）
Sg=
.1O93-I-39.68*k3—26.392-1-6.OilO*K<2.OOOOOOe-OOIj.3^OCQDOOe-001）
3=
2.476-4?
.93*h3+S2.4.6"2-17.til*K<3.OOUDOOe-OO1,4..ODOOOUe-UO1J
xiS
1.30aO0Oe-OOlF5.S31403e-0Ol
三样条第二榮边界峯f牛插直半吉果
S=
.5Q00-I-.380丘764*3：
'3（O,1.OOOOOOs-OOlJ
SS=
5108-9,O7O+k*3+3.247*sc2-t-.5S73e-l*x<1,OOOOOOe-OOU2,OOOOOOe-OO1>
2=
.I6fi3-»-r33,93*x"3-22.55*^c"2^-5.215j»：
x<2-□□□□口口点一（JEJ1,3.0000口□亡一口门1、
S3=
1.SO7—2S.80*1="3-1-32.LI*19-*k（3-000000&—OO1j4-OOOOOOe—OD1>
S
1.SOOOOOe-OOl,,O
5.530212—口Olj.>>
Xi=0.36时
拉搭翎曰栖值绪弟
T.1693S3e-OO1
三样毎第一粪边界号件插値细果
.5Q□□-«-*dODO*K—.3E9B*3E"Z-I-3.09S*h□（Dj1*OOOOOOe—OO1>
S（k>=
“5147—10.99*3e"3-Hlx013*21"2—.4O7Se-1（1*OOOOODe—OO1,2-OOODOOe—OO1>
3（x>=
.1093-1-39.60*m：
*3-26.39+x2+6.O40*x<2.0OOOOOer-OO1,3.OOOOOOer-OO1）
2.475-4T.933+52.46*x2-17.61*k<3,□OOOOOe-OO1a4,OOOOOOt-OOI）
mS
号”0OQQOOe-OQt6t90-sl23O^-QQI
三彳羊朵弟二粪边界祭俏插T亘细乐
-巧QOCH■”38n5*zH-1.75^*x"3“5103—9.OZO+ic"3-1-3x24子狀託"2+_5575e-1（1.□OOOOOe—OOlj2.OOOOOOe—OO1>
“1663-1-33.S3*ie"3—22.SS*3t"2-1-5.21<2.0□OOOOe-O□13.OOOOOOe-OO1）
1.807-26.84*x^3-1-32.14*x*2-l1.19*x<3.000000e-001p4.000000e-001）
x±S
3,600000e-001,O
U.LI
6.315-45Be-OOI^>>|
2将多种插值函数即原函数图像画在同一张图上
0.86
原國数
拉格朗日摘值
牛顿插值摘值
+三样条插值
4.5.1增加插值节点观察误差变化
Jan宜卞熄致为w时旳宜孜果
插值节点数为20时的插值效果
从上面三张图可以看出增加插值节点并不能改善差之效果
4.5.2车门曲线
升级会员 