【正文】 r j=2:n for i=j:n A(i,j)=(A(i,j1) A(i1,j1))/(X(i)X(ij+1))。 end b=poly(X(j1))。q1=conv(q,b)。 c1=c1*j。 q=q1。 end C=A(n,n)。 b=poly(X(n))。 q1=conv(q1,b)。 for k=(n1):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。 return 例 *21給出節(jié)點(diǎn)數(shù)據(jù) f()=， f()= ，f()=， f()= ， f()=， f()=五階牛頓插值多項式和差商，并寫出估計誤差的公式。 解： 輸入程序 X=[ ]。 Y=[ ]。 [A,C,L,wcgs,Cw]= newpoly (X,Y) 運(yùn)行后輸出差商矩陣 A,五階牛頓插值多項式 L及其系數(shù)向量C, 插值余項公式 L及其向量 Cw如下 A = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C = L =7813219284746629/36028797018963968*x^5+583849564517807/9007199254740992*x^4+593245028711263/281474976710656*x^3+3823593773002357/1125899906842624*x^2321902673270315/70368744177664*x+308328649211299/281474976710656 wcgs = 1/720*M*(x^679/25*x^514201/2500*x^4+4934097026900981/281474976710656*x^3+154500712237335/35184372088832*x^28170642380559269/562949953421312*x+5212760744134241/36028797018963968) Cw = 即 L =*x+*x^2+*x^3+*x^*x^5. 估計其誤差的公式為 ? ? ? ? ? ? ? ?? ?? ?? ?? ?? ?!665 ??= xxxxxxfxR ?例 *22求函數(shù) 在 [2,6]上五階牛頓插值多項式，估計其誤差的公式和誤差限公式，用它們計算 f()，并估計其誤差。 ? ? 57 xexf =解： 輸入程序 X=2:4/5:6。 Y=7*exp(X/5)。 [A,C,L,wcgs,Cw]= newpoly(X,Y), x1=2::6。 M=max(7*exp(x1/5)/(5^6)), 運(yùn)行后輸出差商矩陣 A, 五階牛頓插值多項式 L及其系數(shù)向量C, 插值余項公式 L及其向量 Cw如下 A = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 輸入 MATLAB程序 syms x wcgs1=1/720*M*(x^624*x^5+1172/5*x^45952/5*x^3+7276634802928539/2199023255552*x^25237461186650519/1099511627776*x+6085939356447121/2199023255552), 運(yùn)行后輸出誤差限公式 wcgs1如下 wcgs1 =7080139490261079/37778931862957161709568*x^6+21240418470783237/4722366482869645213696*x^52074480870646496147/47223664828696452136960*x^4+658452972594280347/2951479051793528258560*x^351519589424422493159618890033581/83076749736557242056487941267521536*x^2+37081955776313991530044840850001/41538374868278621028243970633760768*x43089299572915358450150147903559/83076749736557242056487941267521536 輸入 MATLAB程序 x=。 y=9721799720875/1152921504606846976*x^53503994098647815/9223372036854775808*x^4+160742021798419/18014398509481984*x^31251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x3940156929554013/562949953421312 wcgs2=1/720*M*(x^624*x^5+1172/5*x^45952/5*x^3+7276634802928539/2199023255552*x^25237461186650519/1099511627776*x+6085939356447121/2199023255552) 運(yùn)行后輸出 的近似值 y，及其誤差限 wcgs2如下 y = wcgs2 = 二、 牛頓插值及其誤差估計的 MATLAB程序 輸入?yún)?shù)： n+1個節(jié)點(diǎn) (xi,yi)(i=1,2,…, n+1)的橫坐標(biāo)向量 X和縱坐標(biāo)向量 Y， x是以向量形式輸入的 m個插值點(diǎn)， M在 [a,b]上滿足 輸出參數(shù)： 向量 x的插值向量 y及其誤差限 R 將程序保存為 M文件，保存到 MATLAB工作目錄下即可使用函數(shù)。 ? ?? ? Mxf n ??1function [y,R]= newcz(X,Y,x,M) n=length(X)。 m=length(x)。 for t=1:m z=x(t)。 A=zeros(n,n)。A(:,1)=Y39。 s=。 p=。 q1=。 c1=。 for j=2:n for i=j:n A(i,j)=(A(i,j1) A(i1,j1))/(X(i)X(ij+1))。 end q1=abs(q1*(zX(j1)))。c1=c1*j。 end C=A(n,n)。q1=abs(q1*(zX(n)))。 for k=(n1):1:1 C=conv(C,poly(X(k)))。d=length(C)。 C(d)=C(d)+A(k,k)。 end y(k)= polyval(C, z)。 end R=M*q1/c1。 return 例 *23已知 sin30 176。 =， sin45 176。 =， sin60 176。 =，用牛頓插值法求 sin40 176。 的近似值，估計其誤差，并與例 *16的計算結(jié)果比較。 解： ? ? ? ? ? ? ? ? ? ? ? ?92s i n40s i n,360,445,630,1c o s,c o s,s i n???? =====?=???=???=???? 求則取 Mxxfxxfxxf方法一（牛頓插值及其誤差估計的 MATLAB主程序） 輸入MATLAB程序 x=2*pi/9。M=1。 X=[pi/6 ,pi/4, pi/3]。 Y=[,]。 [y,R]= newcz(X,Y,x,M) 運(yùn)行后輸出 y = R = 方法二（求牛頓插值多項式和差商的 MATLAB主程序） 輸入 MATLAB程序 x=2*pi/9。 X=[pi/6 ,pi/4, pi/3]。 Y=[,]。 M=1。 [A,C,L,wcgs,Cw]= newpoly(X,Y), y=polyval(C,x) 運(yùn)行后輸出結(jié)果 A = 0 0 0 C = L =1583578379808357/4503599627370496*x^2+1408883391907005/1125899906842624*x132405829044691/2251799813685248 wcgs =1/6*M*(x^33/4*x^2*pi+4012734077357799/2251799813685248*x7757769783530263/18014398509481984) Cw = y = 上述兩種方法計算 y的結(jié)果相同 ,同時與例 *16的結(jié)果也相同。 三、 牛頓插值法的 MATLAB綜合程序 輸入?yún)?shù)： n+1個節(jié)點(diǎn) (xi,yi)(i=1,2,…, n+1)的橫坐標(biāo)向量 X和縱坐標(biāo)向量 Y， x是以向量形式輸入的 m個插值點(diǎn)， M在 [a,b]上滿足 輸出參數(shù)： 向量 x的插值向量 y及其誤差限 R， n次牛頓插值多項式 L及其系數(shù)向量 C， 差商的矩陣 A 將程序保存為 M文件，保存到 MATLAB工作目錄下即可使用函數(shù)。 ? ?? ? Mxf n ??1function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X)。 m=length(x)。 for t=1:m z=x(t)。 A=zeros(n,n)。A(:,1)=Y39。 s=。 p=。 q1=。 c1=。 for j=2:n for i=j:n A(i,j)=(A(i,j1) A(i1,j1))/(X(i)X(ij+1))。 end q1=abs(q1*(zX(j1)))。c1=c1*j。 end C=A(n,n)。q1=abs(q1*(zX(n)))。 for k=(n1):1:1 C=conv(C,poly(X(k)))。 d=length(C)。C(d)=C(d)+A(k,k)。 end y(k)= polyval(C, z)。 end R=M*q1/c1。L(k,:)=poly2sym(C)。 return 例 *24給出節(jié)點(diǎn)數(shù)據(jù) f()=， )= ，f()=， f()= 作三階牛頓插值多項式，計算 f()，并估計誤差。 解： 輸入程序 syms M,X=[4,0,1,2]。 Y =[27,1,2,17]。 x=。 [y,R,A,C,P]=newdscg(X,Y,x,M) 運(yùn)行后輸出插值 y 及其誤差限公式 R，三階牛頓插值多項式 P及其系數(shù)向量 C，差商的矩陣 A如下 y = R = 1323077530165133/562949953421312*M（即 R =*M） A= 0 0 0 0 0 0 C = P = 11/12*x^3+17/4*x^225/6*x+1 例 *25求將區(qū)間 [0,π/2]分成 n等份（ n=2,3），用 y=f(x)=sin(x)產(chǎn)生 n+