龙格函数.doc

对龙格函数在区间-1,1上取n=10,20等距节点，分别作多项式插值、三次样条插值比较各结果。1、 多项式插值：在区间-1,1上取n=10,20等距节点，带入拉格朗日插值多项式中，求出各个节点的插值，并利用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=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end 建立一个polynomial.m文件，用于多项式插值的实现，代码如下：%lagrange插值n=10x=-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)取n=20x=-1:0.1: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)2三次样条插值：所谓三次样条插值多项式是一种分段函数，它在节点分成的每个小区间上是3次多项式，其在此区间上的表达式如下：因此，只要确定了的值，就确定了整个表达式，的计算方法如下：令：则满足如下个方程：对于第一种边界条件下有如果令那么解就可以为 求函数的二阶导数： syms x f=sym(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:1 y=1/(1+25*x2)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.mn=10x=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)得到如下图像：n=20x=linspace(-1,1,21);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)得到如下图像附录三次样条插值主要代码：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-1 u(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-1 A(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:kmax if 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