数值分析作业三次样条插值.docx
- 文档编号:9366821
- 上传时间:2023-02-04
- 格式:DOCX
- 页数:22
- 大小:48.15KB
数值分析作业三次样条插值.docx
《数值分析作业三次样条插值.docx》由会员分享,可在线阅读,更多相关《数值分析作业三次样条插值.docx(22页珍藏版)》请在冰豆网上搜索。
数值分析作业三次样条插值
数值计算方法作业
实验名称
实验4.3三次样条插值函数
(P126)
4.5三次样条插值函数的收敛性
(P127)
实验时间
姓名
班级
学号
成绩
实验4.3三次样条差值函数
实验目的:
掌握三次样条插值函数的三弯矩方法
实验函数:
x
0.0
0.1
F(x)
0.5000
0.5398
0.2
0.3
0.4
0.5793
「0.6179
0.7554
求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)•…
f[X°,Xi,X2,...Xn](X-X°)(X-Xi)...(X
-XnJ)
f[Xi,Xj]
f(Xi)-f(Xj)
Xi-Xj
f[Xi,Xj,Xk]=
其中*.
f[Xj,Xk]-f[K,Xj]
Xk-Xi
3
hi
hiMi4hiMiyi
-6)(6*,皿"]
S(x)二M
i4
(Xi-x)
6hi
式中Mi=S(Xi).
因此,只要确定了Mi的值,就确定了整个表达式,
Mi的计算方法如下:
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 h hihii hi1 hihi1 6 hihii (yii-yi hi屮 上4)=6f[XijXi,Xii]hi 则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 对于第二类边界条件有 ;2M0+人0M1=d0 »nMnJ*2Mn=dn d° 其中 dn 6'0 2(f[x°,X1]-m。 )2(1-0)M0 h1 6"n (m° hn f[Xn4,Xn])2(1-Un)Mn 那么解就可以为 '2.7^0..... 田.2.21 Ym° M1 M2 d0 d1 d2 n4..2..n』 Hn Jn..2 丿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' ); 'r+' plot(a,S, holdon; 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.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 原國数 拉格朗日摘值 牛顿插值摘值 +三样条插值 010.15020.250.303604 0.05 4.5.1增加插值节点观察误差变化 Jan宜卞熄致为w时旳宜孜果 插值节点数为20时的插值效果 从上面三张图可以看出增加插值节点并不能改善差之效果 4.5.2车门曲线
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 作业 三次 样条插值