1设计目的、要求对龙格函数22511)(xx f +=在区间[-1,1]上取10=n 的等距节点，分别作多项式插值、三次样条插值和三次曲线拟合，画出)(x f 及各逼近函数的图形，比较各结果。

2设计原理（1） 多项式插值：利用拉格朗日多项式插值的方法，其主要原理是拉格朗日多项式，即：01,,...,n x x x 表示待插值函数的1n +个节点，()()nn j k k j j k L x y l x y ===∑，其中0,1,...,j n =；011011()...()()...()()()...()...()...()k k n k k k k k k k n x x x x x x x x l x x x x x x x x x -+-+----=----（2） 三次样条插值：三次样条插值有三种方法，在本例中，我们选择第一边界条件下的样条插值，即两端一阶导数已知的插值方法：00'()'S x f = '()'n n S x f =（3）三次曲线拟合：本题中采用最小二乘法的三次多项式拟合。

最小二乘拟合是利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。

在本题中，n= 10，故有11个点，以这11个点的x 和y 值为已知数据，进行三次多项式拟合，设该多项式为23432xi i i i p a a x a x ax =+++，该拟合曲线只需2[]xi i p y -∑的值最小即可。

3采用软件、设备计算机、matlab 软件4设计内容1、多项式插值：在区间[]1,1-上取10=n 的等距节点，带入拉格朗日插值多项式中，求出各个节点的插值，并利用matlab 软件建立m 函数，画出其图形。

在matlab 中建立一个lagrange.m 文件，里面代码如下： %lagrange 函数function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=kp=p*(z-x0(j))/(x0(k)-x0(j)); end ends=p*y0(k)+s; end y(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、三次样条插值：所谓三次样条插值多项式()n S x 是一种分段函数，它在节点i x 011()n n a x x x x b -=<<⋅⋅⋅<<=分成的每个小区间1[,]i i x x -上是3次多项式，其在此区间上的表达式如下：22331111111()[()()]()()666[,]1,2,,.i i i i i i i i i i i i i i ii i h x x h x x S x x x M x x M y M y M h h h x x x i n --------=-+-+-+-∈=⋅⋅⋅，， 因此，只要确定了i M 的值，就确定了整个表达式，i M 的计算方法如下： 令：11111111116()6(,,)i i i i i i i i i i i i i ii i i i i i i h h h h h h y y y y d f x x x h h h h μλμ++++--+++⎧===-⎪++⎪⎨--⎪=-=⎪+⎩，则i M 满足如下1n -个方程：1121,2,,1i i i i i i M M M d i n μλ-+++==⋅⋅⋅-，对于第一种边界条件下有⎪⎪⎩⎪⎪⎨⎧-=+-=+---000110111)'],([62]),['(62h f x x M M h x x f f M M n n n n n n如果令100100016([,]')6('[,])1,,1,,n n n n n n f x x f f f x x d d h h λμ----====那么解就可以为⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----n n n n n n n d d d d M M M M 110110111102222μλμλμλ求函数的二阶导数： >> syms x>> f=sym(1/(1+25*x^2)) f =1/(1+25*x^2) >> diff(f) ans =-(50*x)/(25*x^2 + 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*x^2)endy =得出x=-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1y=0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385编写m文件Three_Spline.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 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)%原曲线 plot(x0,y0,'-r') hold on %三次拟合曲线 plot(x1,y1,'-b')上图中，蓝色部分为三次拟合曲线，红色部分为原曲线6结果分析拉格朗日插值的优点是对于某一区域，不限于被估计点周围，公式简单易实施。

一般认为n 的次数越高，逼近)(x f 的精度就越好，但在本题中，对龙格函数22511)(xx f +=，中间部分插值效果比较好，而对于两端，插值结果是非常不理想的，即龙格现象。

样条函数可以给出光滑的插值曲线，从本题中就能体现出来。

从以上图形可以看出，三次样条插值的图形是比较逼近于原图的，收敛性相对而言是非常好的，但在本题中，仅将原区间分成10个等距区间，因此，逼近效果还不是特别理想，当我们将n 增大时，插值后的曲线越逼近于原曲线。

总的来说，三次样条插值的稳定性比较好，收敛性比较强。

在这三种方法中，三次曲线拟合的效果是最差的，所得的图形与原曲线差距甚远。

最小二乘法中，并不要求拟合后的曲线经过所有已知的点，只需要拟合多项式上的点在某种标准上与定点之间的差距最小即可，因此与原曲线的逼近程度是最差的。

最小二乘法的多项式拟合只适用于多项式，而本题中的函数并不是一个多项式，因此，不建议使用最小二乘法拟合。

参考文献:[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);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=eye(n,n)*2;for k=1:n-1A(k,k+1)=r(k);A(k+1,k)=u(k);endm=A\d;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)&xt(k)<=x(i+1)fx(k)=subs(sx(i),xt(k));endendendplot(xt,fx,'r'); xlabel('x');ylabel('y'); title('f');text(x(fix(n/2)),y(fix(n/2)),'f')hold onplot(x,y,'*')hold off。