插值法-第二次程序題

1、插值法題目1：對Runge函數(shù)在區(qū)間-1,1作下列插值逼近，并和R(x)的圖像進行比較，并對結果進行分析。(1)用等距節(jié)點繪出它的20次Newton插值多項式的圖像。(2)用節(jié)點，繪出它的20次Lagrange插值多項式的圖像。(3)用等距節(jié)點繪出它的分段線性插值函數(shù)的圖像。(4)用等距節(jié)點繪出它的三次自然樣條插值函數(shù)的圖像。程序及分析：(1) 用等距節(jié)點繪出它的20次Newton插值多項式的圖像。Matlab程序如下：%計算均差x=-1:0.1:1;n=length(x);syms zfor i=1:n y(i)=1/(1+25*x(i)*x(i);endN=zeros(n,n); N(:,

2、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:nc(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);endfx=subs(f,z,a);subplot(2,1,1);plot(a,b,'k'

3、,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);endfx=subs(f,z,c);subplot(2,1,2);plot(c,d,'k',c,fx,'r');結果與分析：由下圖可以看出，在區(qū)間-0.6,0.6上，插值多項式可以很好的逼近被插值函數(shù)。而在邊界附近，插值多項式與被插值函數(shù)的差別很大。即出現(xiàn)了Runge現(xiàn)象。主要原因是被插值函數(shù)的任意階導數(shù)不能達到一致有界。其插值余項不趨近零。插值多項式不能收斂到被插值函數(shù)。Runge函數(shù)插值多項式(2)

4、 用節(jié)點，繪出它的20次Lagrange插值多項式的圖像。Matlab程序如下：clear;%插值點for i=1:21 x(i)=cos(2*(i-1)+1)*pi/42);endn=length(x);for i=1:n y(i)=1/(1+25*x(i)*x(i);end%構造插值基函數(shù)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

5、 i=1:n b(i)=1/(1+25*a(i)*a(i);endfx=subs(L,z,a);subplot(2,1,1);plot(a,b,'k',a,fx,'x r');結果與分析：Runge函數(shù)XL插值多項式Runge函數(shù)XL插值多項式Newton插值多項式如下圖所示，使用Chebyshev多項式零點構造的Lagrange插值多項式比較接近原函數(shù)，沒有出現(xiàn)Runge現(xiàn)象。主要原因是其多項式誤差為。(3)用等距節(jié)點繪出它的分段線性插值函數(shù)的圖像。Matlab程序如下：clc;clear;x=-1:0.1:1;n=length(x);syms zfor i=

6、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 onend結果與分析：如下圖所示，分段線性插值多項式比較接近原函數(shù)，沒有出現(xiàn)Runge現(xiàn)象。.分段線性Runge函數(shù)利

7、用線性插值多項式的誤差估計：(4)用等距節(jié)點繪出它的三次自然樣條插值函數(shù)的圖像。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);endfor i=1:n-1h(i)=x(i+1)-x(i);endfor i=1:n-2u(i)=h(i)/(h(i+1)+h(i);r(i)=1-u(i);endG=zeros(n,n);for i=1:n G(i,i)=2;endfor i=2:n-1 G(i,i-1)=u(i-1); G(i,i+1)=r(i-1);endG(n,n-1)=1

8、;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);endsyms 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

9、/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;endfor i=1:n-1 a=x(i):0.01:x(i+1); f=subs(s(i),z,a); plot(a,f,'x r') hold onend結果與分析：三次樣條插值函數(shù)得到的圖像如下：可以看出，三次樣條插值函數(shù)的曲線及其光滑。得到的函數(shù)十分接近被插值函數(shù)。XXXXRunge函數(shù)三次樣條題目2：對函數(shù)：在區(qū)間-1,1作下列插值逼近，并和被插值函數(shù)的圖像進行比較，并對結果進行分析。(1)用等距節(jié)點繪出它的20次Newt

10、on插值多項式的圖像。(2)用節(jié)點，繪出它的20次Lagrange插值多項式的圖像。(3)用等距節(jié)點繪出它的分段線性插值函數(shù)的圖像。(4)用等距節(jié)點繪出它的三次自然樣條插值函數(shù)的圖像。程序及分析：(1)用等距節(jié)點繪出它的20次Newton插值多項式的圖像。Matlab程序如下：clc;clear;%計算均差x=-1:0.1:1;n=length(x);syms z;y=zeros(1,n)for i=1:10 y(i)=sin(pi*x(i);endfor i=11:15 y(i)=cos(pi*x(i);endfor i=15:n y(i)=0;endN=zeros(n,n); N(:,1)

11、=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:nc(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%作圖v=linspace(-1,0,50);u=sin(pi*v);plot(v,u,'k')hold onv=linspace(0,0.5,25);u=cos(pi*v);plot(v,u,'k'

12、;)hold onv=linspace(0.5,1,10000);u=0;plot(v,u,'k')hold ona=-1:0.001:1;fx=subs(f,z,a);plot(a,fx,'r');結果與分析：見下圖被插值函數(shù)插值多項式等距節(jié)點20次Newton插值得到的函數(shù)圖像如下：可以看出，在整個區(qū)間上，插值多項式精度都不是很高。出現(xiàn)了Runge現(xiàn)象。被插值函數(shù)插值多項式(2)用節(jié)點，繪出它的20次Lagrange插值多項式的圖像。Matlab程序如下：clc;clear;%求插值節(jié)點for i=1:21 x(i)=cos(2*(i-1)+1)*pi/42

13、);endn=length(x);y=zeros(1,n);for i=1:n if x(i)<0 y(i)=sin(pi*x(i); elseif x(i)>0.5 y(i)=0; else y(i)=cos(pi*x(i); endend%插值基函數(shù)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;fx=subs(L,z,a);plot(a,fx,'

14、;x r');結果與分析：被插值函數(shù)插值多項式如下圖所示，使用Chebyshev多項式零點構造的Lagrange插值多項式比Newton插值多項式接近原函數(shù)，沒有出現(xiàn)Runge現(xiàn)象。(3)用等距節(jié)點繪出它的分段線性插值函數(shù)的圖像。Matlab程序如下：clc;clear;x=-1:0.1:1;n=length(x);syms z;for i=1:10 y(i)=sin(pi*x(i);endfor i=11:15 y(i)=cos(pi*x(i);endfor i=15:n y(i)=0;end%構造插值多項式for i=1:n-1 l(i)=(z-x(i+1)/(x(i)-x(i+1

15、)*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,'x r') hold onend結果與分析：如下圖所示，分段線性插值多項式比較接近原函數(shù)，沒有出現(xiàn)Runge現(xiàn)象。但是在間斷點處及導數(shù)不存在的點誤差較大。主要是因為這些地方構造的線性函數(shù)斜率較大，不能較好的趨近原函數(shù)。被插值函數(shù)插值多項式(4)用等距節(jié)點繪出它的三次自然樣條插值函數(shù)的圖像。

16、Matlab程序如下：clc;clear;x=-1:0.1:1;n=length(x);syms zfor i=1:10 y(i)=sin(pi*x(i);endfor i=11:15 y(i)=cos(pi*x(i);endfor i=15:n y(i)=0;end for i=1:n-1h(i)=x(i+1)-x(i);endfor i=1:n-2u(i)=h(i)/(h(i+1)+h(i);r(i)=1-u(i);endG=zeros(n,n);for i=1:n G(i,i)=2;endfor i=2:n-1 G(i,i-1)=u(i-1); G(i,i+1)=r(i-1);endG(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);endsyms u v;u=diff(sin(pi*v),v);a=subs(u,v,x(1);b=0;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=