matlab程序等距节点插值多项式,插值法-第二次程序题

插值法-第二次程序题

插值法 题目1：对Runge函数在区间[-1,1]作下列插值逼近，并和R(x)的图像进行比较，并对结果进行分析。 (1)用等距节点绘出它的20次Newton插值多项式的图像。 (2)用节点，绘出它的20次Lagrange插值多项式的图像。 (3)用等距节点绘出它的分段线性插值函数的图像。 (4)用等距节点绘出它的三次自然样条插值函数的图像。 程序及分析： (1) 用等距节点绘出它的20次Newton插值多项式的图像。 Matlab程序如下： %计算均差 x=[-1:0.1:1]; n=length(x); syms z for i=1:n y(i)=1/(1+25*x(i)*x(i)); end N=zeros(n,n); N(:,1)=y ; for j=2:n for k=j:n N(k,j)=(N(k,j-1)-N(k-1,j-1))/(x(k)-x(k-j+1)); end end for t=1:n c(t)=N(t,t) end %构造插值多项式 f=N(1,1); for k=2:n a=1; for r=1:(k-1) a=a*(z-x(r)); end f=f+N(k,k)*a; end %作图 a=[-1:0.001:1]; n=length(a); for i=1:n b(i)=1/(1+25*a(i)*a(i)); end fx=subs(f,z,a); subplot(2,1,1); plot(a,b, k ,a,fx, r ); c=[-0.6:0.001:0.6]; n=length(c); for i=1:n d(i)=1/(1+25*c(i)*c(i)); end fx=subs(f,z,c); subplot(2,1,2); plot(c,d, k ,c,fx, r ); 结果与分析： 由下图可以看出，在区间[-0.6,0.6]上，插值多项式可以很好的逼近被插值函数。而在边界附近，插值多项式与被插值函数的差别很大。即出现了Runge现象。 主要原因是被插值函数的任意阶导数不能达到一致有界。其插值余项不趋近零。插值多项式不能收敛到被插值函数。 Runge函数 插值多项式 (2) 用节点，绘出它的20次Lagrange插值多项式的图像。 Matlab程序如下： clear; %插值点 for i=1:21 x(i)=cos((2*(i-1)+1)*pi/42); end n=length(x); for i=1:n y(i)=1/(1+25*x(i)*x(i)); end %构造插值基函数 syms z; temp=1; for i=1:n lx=1; for j=1:n if i~=j temp=(z-x(j))/(x(i)-x(j)); lx=lx*temp; end end l(i)=lx; end %插值多项式 l=l ; L=y*l; %作图 a=[-1:0.01:1]; n=length(a); for i=1:n b(i)=1/(1+25*a(i)*a(i)); end fx=subs(L,z,a); subplot(2,1,1); plot(a,b, k ,a,fx, x r ); 结果与分析： Runge函数 X L插值多项式 Runge函数 X L插值多项式 Newton插值多项式 如下图所示，使用Chebyshev多项式零点构造的Lagrange插值多项式比较接近原函数，没有出现Runge现象。 主要原因是其多项式误差为。 (3)用等距节点绘出它的分段线性插值函数的图像。 Matlab程序如下： clc;clear; x=[-1:0.1:1]; n=length(x); syms z for i=1:n y(i)=1/(1+25*x(i)*x(i)); end %构造分段线性插值多项式 for i=1:n-1 l(i)=(z-x(i+1))/(x(i)-x(i+1))*y(i)+(z-x(i))/(x(i+1)-x(i))*y(i+1) % l(i)=y(i)+(y(i+1)-y(i))/(x(i+1)-x(i))*(z-x(i)) end %作图 for i=1:n-1 a=[x(i):0.01:x(i+1)]; f=subs(l(i),z,a) plot(a,f, k ) hold on end 结果与分析： 如下图所示，分段线性插值多项式比较接近原函数，没有出现Runge现象。 …… 分段线性 Runge函数 利用线性插值多项式的误差估计： (4)用等距节点绘出它的三次自然样条插值函数的图像。 Matlab程序如下： clc;clear; x=[-1:0.1:1]; n=length(x); syms z; for i=1:n y(i)=1/(1+25*x(i)*x(i)); end for i=1:n-1 h(i)=x(i+1)-x(i); end for i=1:n-2 u(i)=h(i)/(h(i+1)+h(i)); r(i)=1-u(i); end G=zeros(n,n); for i=1:n G(i,i)=2; end for i=2:n-1 G(i,i-1)=u(i-1); G(i,i+1)=r(i-1); end G(n,n-1)=1; G(1,2)=1; d=zeros(1,n); for i=2:n-1 d(i)=6*((y(i+1)-y(i))/h(i)-(y(i)-y(i-1))/h(i-1))/(h(i)+h(i-1)); end syms u v; u=diff(1/(1+25*v*v),v); a=subs(u,v,x(1)); b=subs(u,v,x(n)); d(1)=((y(2)-y(1))/h(1)-a)/h(1)*6; d(n)=(b-(y(n)-y(n-1))/h(n-1))/h(n-1)*6; d=d ; M=inv(G)*d; for i=1:n-1 s(i)=M(i)*(x(i+1)-z)^3/0.6+M(i+1)*(z-x(i))^3/0.6+(y(i)-M(i)*0.01/6)*(x(i+1)-z)/0.1+(y(i+1)-M(i+1)*0.01/6)*(z-x(i))/0.1; end for