三次样条插值.docx
- 文档编号:18418246
- 上传时间:2023-08-16
- 格式:DOCX
- 页数:15
- 大小:107.80KB
三次样条插值.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);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 三次 样条插值