龙格函数

对龙格函数在区间[-1,1]上取n=10,20等距节点，分别作多项式插值、三次样条插值比较各结果。1、多项式插值：在区间[-1,1]上取n=10,20等距节点，带入拉格朗日插值多项式中，求出各个节点的插值，并利用matlab软件建立m函数，画出其图形。在matlab中建立一个lagrange.m文件，里面代码如下：%lagrange函数functiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(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')%插值曲线holdon%原曲线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')%插值曲线holdon%原曲线plot(x0,y1,'-b')-1-0.8-0.6-0.4-0.200.20.40.60.81-60-50-40-30-20-100102三次样条插值：所谓三次样条插值多项式()nSx是一种分段函数，它在节点ix011()nnaxxxxb分成的每个小区间1[,]iixx上是3次多项式，其在此区间上的表达式如下：22331111111()[()()]()()666[,]1,2,,.iiiiiiiiiiiiiiiiihxxhxxSxxxMxxMyMyMhhhxxxin，，因此，只要确定了iM的值，就确定了整个表达式，iM的计算方法如下：令：11111111116()6(,,)iiiiiiiiiiiiiiiiiiiiihhhhhhyyyydfxxxhhhh，则iM满足如下1n个方程：1121,2,,1iiiiiiMMMdin，对于第一种边界条件下有000110111)'],([62]),['(62hfxxMMhxxffMMnnnnnn如果令100100016([,]')6('[,])1,,1,,nnnnnnfxxfffxxddhh那么解就可以为nnnnnnnddddMMMM110110111102222求函数的二阶导数：symsxf=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文件代码如下：forx=-1:0.2:1y=1/(1+25*x^2)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_Spline.mn=10x=linspace(-1,1,11);y=1./(1+25*x.^2);[m,p]=scyt1(x,y,0.0740,-0.0740);holdonx0=-1:0.01:1;y0=1./(1+25*x0.^2);plot(x0,y0,'--b')得到如下图像：-1-0.8-0.6-0.4-0.200.20.40.60.8100.10.20.30.40.50.60.70.80.91xyffn=20x=linspace(-1,1,21);y=1./(1+25*x.^2);[m,p]=scyt1(x,y,0.0740,-0.0740);holdonx0=-1:0.01:1;y0=1./(1+25*x0.^2);plot(x0,y0,'--b')得到如下图像-1-0.8-0.6-0.4-0.200.20.40.60.8100.10.20.30.40.50.60.70.80.91xyfff附录三次样条插值主要代码：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));fork=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;fork=1:n-1A(k,k+1)=r(k);A(k+1,k)=u(k);endm=A\d;ft=d(1);symstfork=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);fori=1:n-1%出点xt对应的y值fork=1:kmaxifx(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')holdonplot(x,y,'*')holdoff