newton_牛顿插值_TheQuestion_源码

clear clc format long n=5; %n=10 %n=20 %n=input('please input the value of n(5,10 or 20):') %第一问 for i=1:n+1 x(i)=-1+2/n*(i-1); y1(i)=f(x(i)); %end x y1 %第二问 %Newton插值多项式 disp('Newton插值多项式'); for j=1:n+1 y(j,1)=f(x(j)); end for j=1:n for i=j:n temp=y(i+1,j)-y(i,j); y(i+1,j+1)=temp/(x(i+1)-x(i+1-j)) ; end c(j)=y(j,j); end c(j+1)=y(j+1,j+1); y %三次样条插值多项式 disp('三次样条插值多项式'); for i=2:n+1 hm(i)=x(i)-x(i-1); end for i=2:n y1(i); u(i-1)=hm(i)/(hm(i)+hm(i+1)); r(i-1)=1-u(i-1); d(i-1)=6/(hm(i)+hm(i+1))*((y1(i+1)-y1(i))/hm(i+1)-(y1(i)-y1(i-1))/hm(i)); end for i=1:n-1 for j=1:n-1 if(i==j) Ms(i,j)=2; elseif(i==j-1) Ms(i,j)=r(i); elseif(i==j+1) Ms(i,j)=u(i); end end end Ms disp('Size of Ms'); size(Ms) %可用追赶法求解 %这里我采用简单一点的方法求得M矩阵 M0=0; M=d*Ms^(-1); %M*Ms %d M(n)=0; hm(1)=(1-(-1))/n; for i=1:n temp=(x(i)-(-1))/(2/n); xs=x(i); if temp<1 S(i)=1/6/hm(1)*((x(2)-xs)^3*M0+(xs-x(1))^3*M(1))+(y1(1)-hm(1)^2*M0/6)*(x(2)-xs)/hm(1)+(y1(2)-hm(1)^2*M(1)/6)*(xs-x(1))/hm(1);; else S(i)=1/6/hm(i)*((x(i+1)-xs)^3*M(i-1)+(xs-x(i))^3*M(i))+(y1(i)-(hm(i))^2*M(i-1)/6)*(x(i+1)-xs)/hm(i)+(y1(i+1)-(hm(i))^2*M(i)/6)*(xs-x(i))/hm(i);; end end S y1 %第三问，直接利用第二问中所求的的差商 %n=5 or 20 for k=1:99 xx(k)=-1+2/100*k; yy(k)=f(xx(k)); end for k=1:99 %xs为一个临时变量 xs=xx(k); for i=1:n+1 if i==1 s(i)=c(i); else s(i)=c(i); for j=1:i-1 s(i)=s(i)*(xs-x(j)); end end end Nn(k)=sum(s); end for i=1:99 temp=(xx(i)-(-1))/(2/n); xs=xx(i); k=ceil(temp); if k==0 k=1; end if k==1 Sn(i)=1/6/hm(1)*((x(2)-xs)^3*M0+(xs-x(1))^3*M(1))+(y1(1)-hm(1)^2*M0/6)*(x(2)-xs)/hm(1)+(y1(2)-hm(1)^2*M(1)/6)*(xs-x(1))/hm(1); else Sn(i)=1/6/hm(k)*((x(k+1)-xs)^3*M(k-1)+(xs-x(k))^3*M(k))+(y1(k)-(hm(k))^2*M(k-1)/6)*(x(k+1)-xs)/hm(k)+(y1(k+1)-(hm(k))^2*M(k)/6)*(xs-x(k))/hm(k); end end xx; yy; Nn; %第4问 ENn=max(abs(yy-Nn)); ESn=max(abs(yy-Sn)); str1='原始曲线,蓝色'; %画图 plot(xx,yy); gtext('原始曲线,蓝色'); hold on plot(xx,Nn,'g'); gtext('牛顿插值，绿色'); plot(xx,Sn,'r') gtext('三次样条插值插值，红色'); hold off