第2章牛顿插值法实现
参考文献：[1]岑宝俊. 牛顿插值法在凸轮曲线修正设计中的应用[J]. 机械工程师,2009,10:54-55.
求牛顿插值多项式和差商的MA TLAB 主程序：
function[A,C,L,wcgs,Cw]=newpoly(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;
（1）保存名为newpoly.m 的M 文件
（2）输入MA TLAB 程序
>> X=[242,243,249,250];
>> Y=[13.681,13.526,13.098,13.095];
>> [A,C,L,wcgs,Cw]=newpoly(X,Y)
输出3阶牛顿插值多项式L 及其系数向量C 差商的矩阵A ，插值余项wcgs 及其
)
()()1(ξ+n n f x R 的系数向量Cw 。

A =
13.6810 0 0 0
13.5260 -0.1550 0 0
13.0980 -0.0713 0.0120 0
13.0950 -0.0030 0.0098 -0.0003
C =
1.0e+003 *
-0.0000 0.0002 -0.0551 4.7634
L =
- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776
wcgs =
(M*(x^4 - 984*x^3 + 363071*x^2 - 59535444*x + 3660673500))/24
Cw =
1.0e+008 *
0.0000 -0.0000 0.0002 -0.0248 1.5253
输入MATLAB程序
>> x=244;
>> y=- (23*x^3)/84000 + (2981*x^2)/14000 - (7757472138947345*x)/140737488355328 + 5237382665812919/1099511627776
y =
13.3976
输入MATLAB程序
>> x=[244,245,246,247,248];
>> y=- (23.*x.^3)./84000 + (2981.*x.^2)./14000 - (7757472138947345.*x)./140737488355328 + 5237382665812919./1099511627776
y =
13.3976 13.2943 13.2143 13.1560 13.1178
第4章 高斯-勒让德积分公式实现
用高斯-勒让德积分公式计算dx e I x ⎰--=112221
π，取代数精度为3和5，再根据截断误差公
式写出误差公式，并将计算结果与精确值进行比较。

（1）建立并保存fun.m 文件命名的M 文件函数
function y=fun(x)
y=exp((-x.^2)./2)./(sqrt(2*pi));
（2）建立并保存GaussR1.m 文件命名的M 文件函数
function[GL,Y,RGn]=GaussR1(fun,X,A)
n=length(X);n2=2*n;Y=feval(fun,X);GL=sum(A.*Y);
sun=1;su2n=1;su2n1=1;wome=1;
syms x
for k=1:n
wome=wome*(x-X(k));
end
wome2=wome^2;Fr=int(wome2,x,-1,1);
for k=1:n2
su2n=su2n*k;
end
syms M
RGn=Fr*M/su2n;
（3）输入程序
代数精度为3：
>> X1=[-1/(3^(1/2)),1/(3^(1/2))];A1=[1,1];
>> [GL1,Y1,Rn1]=GaussR1(@fun,X1,A1)
GL1 =
0.6754
Y1 =
0.3377 0.3377
Rn1 =
M/135
代数精度为5：
>> X2=[-(3/5)^(1/2),0,(3/5)^(1/2)];A2=[5/9,8/9,5/9];
>> [GL2,Y1,Rn2]=GaussR1(@fun,X2,A2)
GL2 =
0.6830
Y1 =
0.2955 0.3989 0.2955
Rn2 =
M/15750
截断误差：
>> syms x
>> fi=int(exp((-x.^2)./2)./(sqrt(2*pi)),x,-1,1);
>> Fs=double(fi),wGL1=double(abs(fi-GL1)),wGL2=double(abs(fi-GL2)) Fs =
0.6827
wGL1 =
0.0073
wGL2 =
3.0777e-004。