实验二:牛顿插值X为各个插值节点,Y为各阶均差
程序如下:
function [A,C,L,wcgs,Cw]= newploy(X,Y)
n=length(X); A=zeros(n,n); A(:,1)=Y';
s=0.0; p=1.0; q=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
b=poly(X(j-1));q1=conv(q,b); c1=c1*j; q=q1;
end
C=A(n,n); b=poly(X(n)); q1=conv(q1,b);
for k=(n-1):-1:1
C=conv(C,poly(X(k))); d=length(C); C(d)=C(d)+A(k,k);
end
L(k,:)=poly2sym(C); Q=poly2sym(q1);
syms M
wcgs=M*Q/c1; Cw=q1/c1;
运行:书P32 例4
X=[0.40 0.55 0.65 0.80 0.90 1.05];
Y=[0.41075 0.57815 0.69675 0.88811 1.02652 1.25382];
[A,C,L,wcgs,Cw]= newploy(X,Y)
输出结果:
A =
Columns 1 through 4
0.4108 0 0 0
0.5782 1.1160 0 0
0.6967 1.1860 0.2800 0
0.8881 1.2757 0.3589 0.1973
1.0265 1.3841 0.4335 0.2130
1.2538 1.5153 0.5249 0.2287
Columns 5 through 6
0 0
0 0
0 0
0 0
0.0312 0
0.0314 0.0003
C =
Columns 1 through 4
0.0003 0.0303 0.1236 0.0296
Columns 5 through 6
0.9901 0.0013
L =
2702819642032443/9223372036854775808*x^5+2181259916290823/72057594037927936*x^4
+4453713600685277/36028797018963968*x^3+2134103104007001/72057594037927936*x^2
+8918190441637495/9007199254740992*x+22964755220417/18014398509481984
wcgs =
1/720*M*(x^6-87/20*x^5+3097/400*x^4-57681/8000*x^3+4166716301494349/1125899906842624* x^2
-4464711044588153/4503599627370496*x+3895001188126157/36028797018963968)
Cw =
Columns 1 through 4
0.0014 -0.0060 0.0108 -0.0100
Columns 5 through 7
0.0051 -0.0014 0.0002