三次样条插值.docx
《三次样条插值.docx》由会员分享，可在线阅读，更多相关《三次样条插值.docx（15页珍藏版）》请在冰豆网上搜索。
三次样条插值
三次样条插值的数值实验
姓名：
王维滨学号：
0842011157姓名：
李佳乐学号：
0842011034
姓名：
谢朝学号：
0842011062姓名：
杨其荣学号：
0842011072
1.实验项目的性质和任务
对三次样条插值进一步理解，并编写matlab程序，实现这些功能
2.算法设计和matlab编程。
总前提：
第i点的横坐标
第i点的纵坐标
的记号
有数值逼近书上的推导，我们令：
，
，
由于未知数的数目多于方程的个数，我们需要增加两个条件才能唯一确定一个分段三次函数
1）D1的三次样条插值
a．实验方案与原理：
我们加上条件：
我们建立三弯矩方程组：
然后采用追赶法迭代求方程组，但是我们在程序中采用简单的方法（矩阵计算）直接求解降低编程难度,
2）D2三次样条插值
3）D3三次样条插值
b．编写程序：
见附录
调用test,每次显示一个图像，关闭后按回车继续显示下一幅。
C．实验结果分析和拓展
（1）
对D1样条：
y11=-1,y1n=1
对D2样条：
对D3样条：
（2）
X=[013.557.510];
Y=[42.51.50.5-0.54];
对D1样条：
y11=2.5,y1n=1.4
对D2样条：
对D3样条：
（3）
（0,1）取五个节点
对D1样条：
y11=0,y1n=0
对D2样条：
对D3样条：
d．实验结论：
（1）f（x）=|x|,
:
D1D2
D3
（2）对于一些型值点
D1D2
:
D3
（3）:
D1:
D2
D3
附录
D1YangTiao.m
functionsvalue=D1YangTiao（xpt,ypt,y11,y1n）
[~,length]=size（xpt）;
u（length）=0;
nemda（length）=0;
u（2:
length-1）=u2_n_1（xpt）;
nemda（2:
length-1）=nemda2_n_1（xpt）;
farr（length）=0;
farr（2:
length-1）=f2_n_1（xpt,ypt）;
farr
（1）=（（ypt
（2）-ypt
（1））/（xpt
（2）-xpt
（1））-y11）/（xpt
（2）-xpt
（1））;
farr（length）=（y1n-（ypt（length）-ypt（length-1））/（xpt（length）-xpt（length-1）））/...
（xpt（length）-xpt（length-1））;
aarr（length,length）=0;
aarr（1,1:
2）=[21];
fori=2:
length-1
aarr（i,i-1:
i+1）=[u（i）,2,nemda（i）];
end
aarr（length,length-1:
length）=[12];
darr（length）=0;
darr=farr*6;
M=solveM（aarr,darr）;
svalue=GetResult（xpt,ypt,M）;
end
D2YangTiao.m
functionsvalue=D2YangTiao（xpt,ypt,y21,y2n）
[~,length]=size（xpt）;
u（length）=0;
nemda（length）=0;
u（2:
length-1）=u2_n_1（xpt）;
nemda（2:
length-1）=nemda2_n_1（xpt）;
farr（length）=0;
farr（2:
length-1）=f2_n_1（xpt,ypt）;
aarr（length-2,length-2）=0;
aarr（1,1:
2）=[2nemda
（2）];
fori=3:
length-2
aarr（i-1,i-2:
i）=[u（i）,2,nemda（i）];
end
aarr（length-2,length-3:
length-2）=[u（length-1）2];
darr（length-2）=0;
darr
（1）=6*farr
（2）-u
（2）*y21;
darr（2:
length-3）=farr（3:
length-2）;
darr（length-2）=farr（length-1）*6-nemda（length-1）*y2n;
M（length）=y2n;
M
（1）=y21;
M（2:
length-1）=solveM（aarr,darr）;
symsx;
svalue（length-1）=x;
fori=1:
length-1
svalue（i）=ypt（i）+（ypt（i+1）-ypt（i））/（xpt（i+1）-xpt（i））*（x-xpt（i））-...
（M（i+1）/6+M（i）/3）*（xpt（i+1）-xpt（i））*（x-xpt（i））+...
M（i）/2*（x-xpt（i））^2+（M（i+1）-M（i））/（xpt（i+1）-xpt（i））/6*（x-xpt（i））^3;
end
D3YangTiao.m
functionsvalue=D3YangTiao（xpt,ypt）
[~,length]=size（xpt）;
xn_1=xpt（length）+xpt
（2）-xpt
（1）;
u（length）=0;
nemda（length）=0;
u（2:
length-1）=u2_n_1（xpt）;
u（length）=（xpt（length）-xpt（length-1））/（xn_1-xpt（length-1））;
nemda（2:
length-1）=nemda2_n_1（xpt）;
nemda（length）=（xn_1-xpt（length））/（xn_1-xpt（length-1））;
farr（length）=0;
farr（2:
length-1）=f2_n_1（xpt,ypt）;
farr（length）=（（ypt
（2）-ypt（length））/（xpt
（2）-xpt
（1））-（ypt（length）-ypt（length-1））/（xpt（length）-xpt（length-1）））/（xn_1-xpt（length-1））;
aarr（length-1,length-1）=0;
aarr（1,1:
2）=[2nemda
（2）];
aarr（1,length-1）=u
（2）;
fori=3:
length-1
aarr（i-1,i-2:
i）=[u（i）,2,nemda（i）];
end
aarr（length-1,length-2:
length-1）=[u（length）2];
aarr（length-1,1）=nemda（length）;
darr（length-1）=0;
darr（1:
length-1）=6*farr（2:
length）;
M（length）=0;
M（2:
length）=solveM（aarr,darr）;
M
（1）=M（length）;
svalue=GetResult（xpt,ypt,M）;
end
f2_n_1.m
functionfarray=f2_n_1（xpt,ypt）
[~,length]=size（xpt）;
fori=2:
length-1
farray（i-1）=（（ypt（i+1）-ypt（i））/（xpt（i+1）-...
xpt（i））-（ypt（i）-ypt（i-1））/（xpt（i）-xpt（i-1）））/...
（xpt（i+1）-xpt（i-1））;
end
end
u2_n_1.m
functionuarr=u2_n_1（xpt）
[~,length]=size（xpt）;
uarr（length-2）=0;
fori=2:
length-1
uarr（i-1）=（xpt（i）-xpt（i-1））/（xpt（i+1）-xpt（i-1））;
end
end
nemda2_n_1.m
functionnemdaarr=nemda2_n_1（xpt）
[~,length]=size（xpt）;
uarr（length-2）=0;
fori=2:
length-1
nemdaarr（i-1）=（xpt（i+1）-xpt（i））/（xpt（i+1）-xpt（i-1））;
end
end
getIndex.m
functionindex=getIndex（xpt,xx）
[~,length]=size（xpt）;
ifxx（1）||xx>xpt（length）
index=-1;
else
fori=1:
length
ifxx<=xpt（i）
index=i-1;
break;
end
end
end
end
PlotYangTiao.m
functionPlotYangTiao（xpt,strFunArr）
[~,length]=size（xpt）;
fori=1:
length-1
xx=xpt（i）:
0.0001:
xpt（i+1）;
yy=subs（strFunArr（i）,xx）;
plot（xx,yy）;holdon
end
end
YangTiao.m
function[result,strFunArr]=YangTiao（xpt,ypt,xx,paraArr）
[~,length]=size（xpt）;
ifxx（1）||xx>xpt（length）
result=0;
strFunArr=0;
disp"±ØÐëÔÚÇø¼äÄÚ"
return;
end
strFunArr（length）=0;
ifparaArr
（1）==1%¼ÆËãD1ÑùÌõ
strFunArr=D1YangTiao（xpt,ypt,paraArr
（2）,paraArr（3））;
elseifparaArr
（1）==2%¼ÆËãD2ÑùÌõ
strFunArr=D2YangTiao（xpt,ypt,paraArr
（2）,paraArr（3））;
elseifparaArr
（1）==3%¼ÆËãD3ÑùÌõ
strFunArr=D3YangTiao（xpt,ypt）;
end
fori=1:
length
ifabs（xx-xpt（i））<=10^（-7）
result=ypt（i）;
return
end
end
index=getIndex（xpt,xx）;
result=subs（strFunArr（index）,xx）;
end
solveM.m
functionM=solveM（aarr,darr）
M=inv（aarr）*darr';
end
comparePlot.m
functioncomparePlot（xvec,yvec,funarr）
PlotYangTiao（xvec,funarr）;holdon
plot（xvec,yvec）;
end
GetResult.m
functionsvalue=GetResult（xpt,ypt,M）
symsx;
[~,length]=size（xpt）;
svalue（length-1）=x;
fori=1:
length-1
svalue（i）=ypt（i）+（ypt（i+1）-ypt（i））/（xpt（i+1）-xpt（i））*（x-xpt（i））-...
（M（i+1）/6+M（i）/3）*（xpt（i+1）-xpt（i））*（x-xpt（i））+...
M（i）/2*（x-xpt（i））^2+（M（i+1）-M（i））/（xpt（i+1）-xpt（i））/6*（x-xpt（i））^3;
end
test.m
%
（1）
xvec=-1:
0.5:
1;
yvec=abs（xvec）;
%¶ÔD1
[~,funarr]=YangTiao（xvec,yvec,0.3,[1-11]）;
comparePlot（xvec,yvec,funarr）;
pause
%¶ÔD2
[~,funarr]=YangTiao（xvec,yvec,0.3,[200]）;
comparePlot（xvec,yvec,funarr）;
pause
%¶ÔD3
[~,funarr]=YangTiao（xvec,yvec,0.3,[3]）;
comparePlot（xvec,yvec,funarr）;
pause
%
（2）
xvec=[013.557.510];
yvec=[42.51.50.5-0.54];
pause
%D1
[~,funarr]=YangTiao（xvec,yvec,3,[12.51.4]）;
comparePlot（xvec,yvec,funarr）;
pause
%D2
[~,funarr]=YangTiao（xvec,yvec,3,[200]）;
comparePlot（xvec,yvec,funarr）;
pause
%D3
[~,funarr]=YangTiao（xvec,yvec,3,[3]）;
comparePlot（xvec,yvec,funarr）;
pause
%（3）
xvec=0.001:
0.2:
1;
yvec=sin（xvec）./xvec;
%D1
[~,funarr]=YangTiao（xvec,yvec,0.3,[100]）;
comparePlot（xvec,yvec,funarr）;
pause
%d2
[~,funarr]=YangTiao（xvec,yvec,0.3,[200]）;
comparePlot（xvec,yvec,funarr）;
pause
%d3
[~,funarr]=YangTiao（xvec,yvec,0.3,[3]）;
comparePlot（xvec,yvec,funarr）;
升级会员 