【正文】 作為誤差估計上限。選擇要計算的 x 所在的區(qū)間的端點，插值效果較好。 li(x) 每個 li 有 n 個根 x0 … xi … xn ? = = = n j j ? i j i n i i i x x C x x x x x x C x l 0 0 ) ( ) )...( )...( ( ) ( ? = = j ? i j i i i i x x C x l ) ( 1 1 ) ( ?=? = njij jiji xxxxxl0)()()(?== niiin yxlxL0)()( 拉格朗日 多項式 與 有關(guān)，而與 無關(guān) 節(jié)點 f 拉格朗日插值公式 Lagrange插值 公式 (利用 插值基函數(shù) 很容易得到 ): 含義直觀 ,結(jié)構(gòu)緊湊 ,在理論分析中非常方便 。我們可以有具有“ 承襲性 ”的所謂牛頓公式。 如輸入程序 Y = poly (3)；則返回 1 3，即其特征多項式為 x3 如 對多項式 ， 計算在 x=5,7,9的值。V39。y=poly2sym(c) ；則返回 y =x3 如輸入程序 c=[1 3]。C=conv(A,B) 則返回 C = 1 5 6 A、 B為兩個多項式的系數(shù) (五 ) DECONV 函數(shù) 調(diào)用格式 ： [Q,R] =deconv (B,A) conv函數(shù)的逆函數(shù)，返回兩個多項式相除的多項式系數(shù)及其余數(shù)，即反卷積和 若輸入程序 C=[1 5 7]。[B,R]=deconv(C,A) 則輸出 B = 1 2 R =0 0 0 (六 ) roots(poly(1:n))命令 調(diào)用格式： roots(poly(1:n)) (七 ) det(a*eye(size (A)) A)命令 調(diào)用格式： b=det(a*eye(size (A)) A) 如輸入程序 C=[1 3 2]。 l01= poly(X(2))/( X(1) X(2)), l11= poly(X(1))/( X(2) X(1)), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=。x=0::1。x=pi/6。x=pi/6。 p1=poly(X(1))。 l01= conv ( conv (p2, p3), p4)/(( X(1) X(2))* ( X(1) X(3)) * ( X(1) X(4))), l11= conv ( conv (p1, p3), p4)/(( X(2) X(1))* ( X(2) X(3)) * ( X(2) X(4))), l21= conv ( conv (p1, p2), p4)/(( X(3) X(1))* ( X(3) X(2)) * ( X(3) X(4))), l31= conv ( conv (p1, p2), p3)/(( X(4) X(1))* ( X(4) X(2)) * ( X(4) X(3))), l0=poly2sym (l01), l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31), P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4), 運行后輸出基函數(shù) l0， l1， l2和 l3及其插值多項式的系數(shù)向量 P為 l0 =1/24*x^3+1/8*x^21/12*x， l1 =1/4*x^31/4*x^2x+1 l2 =1/3*x^3+4/3*x， l3 =1/8*x^3+1/8*x^21/4*x P = 1 4 4 1 輸入程序 L=poly2sym (P),x=。 L=ones(m,m)。 l(k,:)=poly2sym (V) end C=Y*L1。 Y=[ ]。s=。 for j=1:n if j~=k p=p*(zX(j))/(X(k)X(j))。 end y(i)=s。 =， sin60 176。 M=1。 )= … ) ( ) ( 0 0 1 0 1 0 1 x x x x y y y x P ? = ) ( 0 0 1 0 1 0 x x x x x x y ? = ( ( ) ) f f f[x0,x1] 二次牛頓插值多項式 我們再看 線性插值 的 點斜式 : ) ( 0 0 x x y ? = f[x0,x1] 常數(shù) (差商 ) 由此啟發(fā)，我們希望二次插值也能類似地有有規(guī)律的組合表達(dá)式 : P2(x)=?0 + ?1(xx0) + ?2(xx0)(xx1) 利用 P2(x0)=y0有 : ?0 = y0 , 利用 P2(x1)=y1有 : ?1 = 0 1 0 1 x x x x ( ( ) ) f f = f[x0,x1] , 利用 P2(x2)=y2有 : ?2 = f[x0,x1] (x2x0)(x2x1) (x2x0)(x2x1) 0 x x2 ( ( ) ) f f (x2x0) f[x0,x2] f[x0,x1] x2 x1 = = f[x0,x1,x2] 。 k 階均差可表示為函數(shù)值 f(x0), f(x1),…, f(xk)的線性組合 , 即 f[x0,x1,…,x k]= f(xj) (xjxj+1)… (xjxk) …( xjxj1) (xjx0) ∑ k j=0 這個性質(zhì)可用歸納法證明 . 這個性質(zhì)也表明均差與節(jié)點的排列次序無關(guān) ,稱為均差的對稱性 ,即 f[x0,x1,…,x k]= f[x1,x0,x2,…,x k]=… = f[x1, …, x k ,x0] f[x0, x1,…,x k] = f(n)(ξ) n! ],[ ba?? 3186。近似計算時，牛頓插值法計算量小，當(dāng)增加節(jié)點時只需增加一項，因而比較方便。本節(jié)給出具有承襲性的牛頓插值法的 MATLAB實現(xiàn)程序。 A=zeros(n,n)。 q=。q1=conv(q,b)。 b=poly(X(n))。 C(d)=C(d)+A(k,k)。 Cw=q1/c1。 [A,C,L,wcgs,Cw]= newpoly (X,Y) 運行后輸出差商矩陣 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()，并估計其誤差。 M=max(7*exp(x1/5)/(5^6)), 運行后輸出差商矩陣 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), 運行后輸出誤差限公式 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=。 for t=1:m z=x(t)。 p=。 end q1=abs(q1*(zX(j1)))。 for k=(n1):1:1 C=conv(C,poly(X(k)))。 end R=M*q1/c1。 =，用牛頓插值法求 sin40 176。 X=[pi/6 ,pi/4, pi/3]。 Y=[,]。 ? ?? ? Mxf n ??1function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X)。A(:,1)=Y39。 c1=。 end C=A(n,n)。C(d)=C(d)+A(k,k)。 return 例 *24給出節(jié)點數(shù)據(jù) f()=， )= ，f()=， f()= 作三階牛頓插值多項式，計算 f()，并估計誤差。 [y,R,A,C,P]=newdscg(X,Y,x,M) 運行后輸出插值 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+1