资源描述:
《数值分析三次样条插值函数》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、数值分析第三次程序作业PB09001057孙琪【问题】对函数fx=ex,x∈[0,1]构造等距节点的三次样条插值函数,对以下两种类型的样条函数1.三次自然样条2.满足S'0=1,S'1=e的样条并计算如下误差:maxi{fxi-12-Sxi-12,i=1,…,N}这里xi-12为每个小区间的中点。对N=10,20,40比较以上两组节点的结果。讨论你的结果。【三次样条插值】在每一个区间[t1,t2],…,[tn-1,tn]上,S都是不同的三次多项式,我们把在[ti-1,ti]上表示S的多项式记为Si,从而,Sx=S0xx∈[t0,t1]S1
2、xx∈[t1,t2]…Sn-1xx∈[tn-1,tn]通过在节点处函数值、一阶导数和二阶导数的连续性可以得到:Si-1ti=yi=Siti1≤i≤n-1Si-1'ti=Si'tilimx→ti+S''x=zi=limx→ti-S''(x)再给定z0和zn的值就构成了4n个条件,而三次样条插值函数共4n个系数,故可以通过这4n个条件求解三次样条函数的系数,从而求得该三次样条插值函数。特别的,当z0=zn=0时称为自然三次样条。一、自然三次样条插值【自然三次样条插值算法】1.由上面的分析可知,求解三次样条函数实际上就是求解一个矩阵:u1h1h
3、1u2h2h2u3h3………hn-3un-2hn-2hn-2un-1z1z2z3…zn-2zn-1=v1v2v3…vn-2vn-1其中hi=ti+1-ti,ui=2(hi+hi-1),ui=6hi(yi+1-yi),vi=bi-bi-1所以自然三层次样条插值的算法就是在得到端点的函数值,一次导数值和二次导数值,然后根据上述求解矩阵得到v,代入自然三次样条的表达式即可。2.根据题目中所给出的误差估计,计算在区间中点处的最大误差。【实验】通过Mathematica编写程序得到如下结果:N=101.计算得到zi的值为:由此可以得到各个区间的自然
4、三次样条插值函数。2.计算得到各区间中点的最大误差为:自然三次样条插值函数与原函数的图像为:可以看出两者基本吻合,与我们计算得到的最大误差符合,实验结果成立。【结果】同样的,我们可以得到该题目的结果:NMaxErrorofgrid(1)100.00124198837692200.000310817169505400.0000777244871042通过观察不难看出:e110≈22×e120≈22×e140【分析】通过结果可以看出,区间越分越细时,区间中点的最大误差越来越小,说明自然三次样条插值越来越逼近原来的函数。一、普通三次样条插值【普
5、通三次样条插值算法】1.在给定了S’(0)=1和S’(1)=e,类似于自然三次样条插值,但此时的矩阵有变化,因此在计算时与之前的自然三次样条插值有所区别。2.根据题目中所给出的误差估计,计算在区间中点处的最大误差。【实验】N=101.计算得到zi的值为:2.计算得到各区间中点的最大误差为:6.95586472865×10-7普通三次样条插值函数与原函数的图像为:可以看出两者几乎吻合,与我们计算得到的最大误差符合,实验结果成立。【结果】同样的,我们可以得到该题目的结果:NMaxErrorofgrid(2)106.95586472865×10
6、-7204.38712914885×10-8402.75377542991×10-9通过观察不难看出:e210≈24×e220≈24×e240【分析】通过结果可以看出,区间越分越细时,区间中点的最大误差越来越小,说明普通三次样条插值越来越逼近原来的函数。【对比】NMaxErrorofgrid(1)MaxErrorofgrid(2)100.001241988376926.95586472865×10-7200.0003108171695054.38712914885×10-8400.00007772448710422.75377542991
7、×10-9通过对比可以看出:同样的三次插值函数,在中点处的误差普通三次插值要比自然三次插值小很多,虽然我们不能确定的说普通三次插值就要比自然三次插值准确很多,但不难猜测到结果跟上述论断差不多。经过对不同函数用两种方法检验,两者的精度并不一定是这种情况,只是对该题有这个结果而已。通过观察不难看出:e110≈22×e120≈22×e140e210≈24×e220≈24×e240【Mathematica程序】自然三次样条插值:f1[x_]:=Exp[x];M=10;t=-1/M;For[i=0,i£M,i++,t=t+1/M;a[i]=t;y[
8、i]=f1[t]];For[i=0,i£M-1,i++,h[i]=a[i+1]-a[i];c[i]=6(y[i+1]-y[i])/h[i]];u[1]=2(h[0]+h[1]);v[1]=c[