Matlab主程序
X=[0.2,0.4,0.6,0.8,1.0];
Y=[0.98,0.92,0.81,0.64,0.38];
M=1;
m=[0,1,11,10];
a=length(m);
x=zeros(size(m));%生成和m相同的零矩阵
for i=1:a
x(i)=0.2+0.08*m(i);
end
x1=linspace(0.2,1.0,10);
[y R A]=Newton(X,Y,x1,M);
[y1 R1 A1]=Newton(X,Y,x,M);
s=interp1(X,Y,x1,'spline');
s1=interp1(X,Y,x,'spline');
plot(X,Y,'or',x1,y,'-b',x1,s,'-k',x,y1,'sb',x,s1,'>k');
legend('Raw','Newton line','Spline line','Newton dots','Spline dots');
Newton插值法函数
function [y,R,A]= Newton(X,Y,x,M)%其中A是如书中31页所示的均差表，理清楚其中的数字关系，就能完成编程
n = length(X); m = length(x);
for t=1:m
z=x(t);
A=zeros(n,n);
A(:,1)=Y';
q1=1.0;c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));
end
q1=abs(q1*(z-X(j-1)));
c1=c1*j;
end
s=A(1,1);
for k=2:n
C=poly(X(1:k-1));
s=s+A(k,k)*polyval(C,z);%虽然看上去思路清晰了，但是感觉不是效率最高的，没有用到比较高级的数学工具。
end
y(t)=s;
R(t)=M*q1/c1;
end
结果图
1.jpg
重点说明
- poly(a1,a2,...,an)表示的是(x-a1)(x-a2)...(x-an)的多项式，从高次到常数项的系数矩阵。
- polyval(A,z)表示的是，以A行列式为系数的多项式，x=z时的取值。
- plot函数中，单引号中的语句可以调整最后呈现的方式。可以通过在主窗口输入help plot进行查阅。
- interp1(X,Y,x,'method')是matlab自带的插值函数。其中X,Y表示插值点，x表示自变量，‘method’共有一下几种插值方式：
'nearest' - nearest neighbor interpolation
'linear' - linear interpolation
'spline' - piecewise cubic spline interpolation (SPLINE)
'pchip' - shape-preserving piecewise cubic interpolation
'cubic' - same as 'pchip'
'v5cubic' - the cubic interpolation from MATLAB 5, which does not
extrapolate and uses 'spline' if X is not equally
spaced.
- linspace(a,b,n)表示的是在[a,b]区间内，插入n个等距节点。
- matlab的数组记数是从1开始的！
