Newton插值和三次样条插值的程序代码.docVIP

Newton插值的程序代码： %求Newton插值多项式 function new=newton(m,y) n=length(m); b=1;new=[zeros(1,n-1),y(1)]; for i=2:n y(i:n)=(y(i:n)-y(i-1:n-1))./(m(i:n)-m(1:n-i+1)); %求差商 b=conv(b,[1,-m(i-1)]); a=[zeros(1,n-length(b)),b]; new=new+y(i)*a; end 三次样条插值的程序代码： %求三次样条插值多项式 function s=myspline(m,y) syms x; n=length(m); M=y; M(2:n)=(M(2:n)-M(1:n-1))./(m(2:n)-m(1:n-1)); M(3:n)=(M(3:n)-M(2:n-1))./(m(3:n)-m(1:n-2)); M=[0,M(3:n),0]; h=m(2:n)-m(1:n-1); c=h(2:n-1)./(h(2:n-1)+h(1:n-2)); a=1-c;b=2+zeros(1,n);M=6*M; c=[0,c];a=[a,0]; M=chase(a,b,c,M); %用追赶法解方程组 for i=1:n-1 %分区间计算插值多项式 s(i)=M(i)*(m(i+1)-x)^3+M(i+1)*(x-m(i))^3+(6*y(i)-M(i)*h(i)^2)*(m(i+1)-x)+(6*y(i+1)-M(i+1)*h(i)^2)*(x-m(i)); s(i)=s(i)/(6*h(i)); end 主程序为： clear;clc; n=5; %指定n的值 syms x %符号数x f=1./(1+25*x^2); %定义函数f(x) i=0:n;m=-1+2*i/n; y=subs(f,m); %计算函数值f(x) disp([When n=,num2str(n),,the values of f(x) at the points of x(i) are:]); disp(y); %求Newton插值多项式 new=newton(m,y);new=round(10000*new)/10000; nout=poly2sym(new);nout=vpa(nout,4); %整理多项式并输出 disp([When n=,num2str(n),,N(x)=]);disp(nout); %求三次样条插值多项式 s=myspline(m,y); for i=1:n %分区间计算插值多项式 a=sym2poly(s(i));b=poly2sym(a);b=vpa(b,4); %整理多项式并输出 disp([In the range of [,num2str(m(i)),,,num2str(m(i+1)),],S(x)=]); disp(b); end k=0:100;xk=-1+0.02*k; %计算xk的值 ll1=subs(f,xk);ll2=polyval(new,xk); %计算f(xk),N(xk) for i=1:length(xk) r=find(m(2:n+1)=xk(i),1,first); ll3(i)=subs(s(r),xk(i)); %计算S(xk) end en=max(abs(ll1-ll2));es=max(abs(ll1-ll3)); %计算E(Nn),E(Sn),并输出 disp([E(Nn)=,num2str(en)]);disp([E(Sn)=,num2str(es)]); plot(xk,ll1,xk,ll2,r,xk,ll3,g);axis([-1,1,-0.1,1.1]); title(Interpolation);legend(Original function,Newton Interpolation,Spline) grid; %将原函数，Newton插值函数，三次样条插值函数画成图形