三次样条拟合范例

1、.1 设计目的 、要求1对龙格函数f ( x)125x 2 在区间 -1,1 上取 n 10 的等距节点 ，分别作多项式插值、三次样条插值和三次曲线拟合，画出 f ( x) 及各逼近函数的图形 ，比较各结果 。2 设计原理（1） 多项式插值 ：利用拉格朗日多项式插值的方法 ，其主要原理是拉格朗日多项式，即：x0 , x1,., xn 表示待插值函数的n1个节点，nLn (x j )yklk ( xj )yj，其中 j0,1,., n ；k0lk (x)( x x0 ).( xxk 1 )( x xk1 ).( x xn )x0 ).(xkxk1 ).( xkxk 1 ).( xk xn )(x

2、k（2） 三次样条插值 ：三次样条插值有三种方法 ，在本例中 ，我们选择第一边界条件下的样条插值 ，即两端一阶导数已知的插值方法 ：S '(x0 )f ' 0S '(xn )f 'n（3）三次曲线拟合 ：本题中采用最小二乘法的三次多项式拟合。最小二乘拟合是利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。 在本题中 ， n=10，故有 11 个点，以这 11 个点的 x和 y 值 为 已 知 数 据 ， 进 行 三 次 多 项 式 拟 合 ， 设 该 多 项 式 为 pxi a4 a3 xi a2 xi 2 axi 3 ，该拟合

3、曲线只需 pxi yi 2 的值最小即可 。3 采用软件 、设备计算机、matlab 软件.专业专注.4 设计内容1、 多项式插值 ：在区间1,1 上取 n10 的等距节点 ，带入拉格朗日插值多项式中，求出各个节点的插值，并利用 matlab 软件建立 m 函数，画出其图形 。在 matlab 中建立一个 lagrange.m 文件，里面代码如下 ：%lagrange函数function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j=kp=p*(z-x

4、0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;end.专业专注.建立一个 polynomial.m文件，用于多项式插值的实现 ，代码如下 ：%lagrange 插值x=-1:0.2:1;y=1./(1+25*x.2);x0=-1:0.02:1;y0=lagrange(x,y,x0);y1=1./(1+25*x0.2);plot(x0,y0,'-r')%插值曲线hold on%原曲线plot(x0,y1,'-b')运行 duoxiangshi.m文件，得到如下图形 ：.专业专注.2、 三次样条插值 ：所谓三次样条插值多

5、项式 Sn ( x) 是一种分段函数 ，它在节点xi (ax0 x1xn 1xn b) 分成的每个小区间 xi 1 , xi 上是 3 次多项式 ，其在此区间上的表达式如下：S( x)13M i 13M i ( yi 1hi2xix( yihi2x xi 1，( xix)( x xi 1)6M i 1)hi6M i )6hihix，, n. xi 1, xi i 1,2,因此，只要确定了 M i 的值，就确定了整个表达式 ， M i的计算方法如下 ：令：.专业专注.ihi， ihi 11 ihi 1hihi 1hidi6( yi 1yiyiyi 1 ) 6 f ( xi 1 , xi , xi

6、 1)hihi 1hi1hi则 M i 满足如下 n1个方程：i M i 12M ii M i 1di，i1,2, n1对于第一种边界条件下有M n 12M n6( f 'n f xn 1 , xn )hn 12M 0M 16( x1 , x0 f ' 0 )h0如果令 0 1,d06( f x1 , x0 f '0 ) ,n 1,dn6( f 'nf xn 1, xn ), 那么解就可以为h0hn 120M 0d 0121M 1d1n 12n 1M n 1d n 1n2M nd n求函数的二阶导数 ：>> syms x>> f=sym(

7、1/(1+25*x2)f =1/(1+25*x2)>> diff(f)ans =-(50*x)/(25*x2 + 1)2将函数的两个端点 ，代入上面的式子中 ：.专业专注.f (-1)= 0.0740f (1)=-0.0740求出从 -1 到 1 的 n=10 的等距节点 ，对应的 x， y 值对应 m 文件代码如下 ：for x=-1:0.2:1y=1/(1+25*x2)endy =得出x=-1-0.8-0.6-0.4-0.200.20.40.60.81y=0.03850.05880.10.20.510.50.20.10.05880.0385编写 m 文件 Three_Splin

8、e.mx=linspace(-1,1,11);y=1./(1+25*x.2);m,p=scyt1(x,y,0.0740,-0.0740);hold onx0=-1:0.01:1;y0=1./(1+25*x0.2);plot(x0,y0,'-b')得到如下图像 ：.专业专注.其中蓝色曲线为原图 ，红色曲线为拟合后的图像。3、 三次曲线拟合 ：这里我们使用最小二乘法的3 次拟合建立一个 Three_fitting .m文件，代码如下 ：%主要代码x=-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1;y=0.0385 0.0588 0.1 0.2

9、0.5 1 0.5 0.2 0.1 0.0588 0.0385;a=polyfit(x,y,3);x1=-1:0.01:1;y1=a(4)+a(3)*x1+a(2)*x1.2+a(1)*x1.3;x0=-1:0.01:1;y0=1./(1+25*x0.2).专业专注.1 25x 21.%原曲线plot(x0,y0,'-r')hold on%三次拟合曲线plot(x1,y1,'-b')上图中，蓝色部分为三次拟合曲线，红色部分为原曲线6 结果分析拉格朗日插值的优点是对于某一区域，不限于被估计点周围 ，公式简单易实施 。一般认为 n 的次数越高 ，逼近 f ( x)

10、的精度就越好 ，但在本题中 ，对龙格函数 f ( x)，中间部分插值效果比较好 ，而对于两端 ，插值结果是非常不理想的 ，即龙格现象 。 样条函数可以给出光滑的插值曲线 ，从本题中就能体现出来 。从以上图形可以看出 ，三次样条插值的图形是比较逼近于原图的 ，收敛性相对而言是非常好的 ，但在本题中 ，仅将原区间分成 10 个等距区间 ，因此，逼近效果还不是特别理想 ，当我们将 n 增大时，插值.专业专注.后的曲线越逼近于原曲线 。总的来说 ，三次样条插值的稳定性比较好 ，收敛性比较强 。在这三种方法中 ，三次曲线拟合的效果是最差的 ，所得的图形与原曲线差距甚远 。最小二乘法中 ，并不要求拟合后的

11、曲线经过所有已知的点 ，只需要拟合多项式上的点在某种标准上与定点之间的差距最小即可 ，因此与原曲线的逼近程度是最差的 。最小二乘法的多项式拟合只适用于多项式 ，而本题中的函数并不是一个多项式 ，因此，不建议使用最小二乘法拟合 。参考文献 :1 李庆扬 王能超等 .数值分析 M. 清华大学出版社2吴振远 .科学计算实验指导书 基于 MATLAB 数值分析 M. 中国地质大学出版社3宋叶志 . MATLAB 数值分析与应用 M. 机械工业出版社 , 2009.07附录三次样条插值主要代码：function m,p=scyt1(x,y,df0,dfn)n=length(x);r=ones(n-1,1

12、);u=ones(n-1,1);d=ones(n,1);r(1)=1;d(1)=6*(y(2)-y(1)/(x(2)-x(1)-df0)/(x(2)-x(1);u(n-1)=1;d(n)=6*(dfn-(y(n)-y(n-1)/(x(n)-x(n-1)/(x(n)-x(n-1);for k=2:n-1u(k-1)=(x(k)-x(k-1)/(x(k+1)-x(k-1); r(k)=(x(k+1)-x(k)/(x(k+1)-x(k-1); d(k)=6*(y(k+1)-y(k)/(x(k+1)-x(k)-(y(k)-y(k-1)/(x(k)-x(k-1)/(x(k+1)-x(k-1);endA=

13、eye(n,n)*2;.专业专注.for k=1:n-1A(k,k+1)=r(k);A(k+1,k)=u(k);endm=Ad;ft=d(1);syms tfor k=1:n-1%求 s(x)即插值多项式p(k,1)=m(k)/(6*(x(k+1)-x(k);p(k,2)=m(k+1)/(6*(x(k+1)-x(k);p(k,3)=(y(k)-m(k)*(x(k+1)-x(k)2/6)/(x(k+1)-x(k);p(k,4)=(y(k+1)-m(k+1)*(x(k+1)-x(k)2/6)/(x(k+1)-x(k);sx(k)=p(k,1)*(x(k+1)-t)3+p(k,2)*(t-x(k)3+p(k,3)*(x(k+1)-t)+p(k,4)*(t-x(k);endkmax=1000;xt=linspace(x(1),x(n),kmax);for i=1:n-1%出点 xt 对应的 y 值for k=1:kmaxif x(i)<=xt(k)&