基于MATLAB軟件的牛頓插值法理論研究

周建杰 禹仁貴

摘 要：牛頓插值法是數(shù)值計算中較為重要的方法之一。本文介紹了牛頓插值法的理論，討論了其優(yōu)缺點及與泰勒展式的關(guān)系。通過一個具體的實例，應(yīng)用MATLAB軟件，直觀地展示了牛頓插值法的實現(xiàn)及其特點。

關(guān)鍵詞：牛頓插值法;泰勒展式;龍格現(xiàn)象;MATLAB

中圖分類號：G642文獻標(biāo)識碼：A文章編號：1003-5168（2018）32-0010-03

Newton Interpolation Method Theory and MATLAB Realization

ZHOU Jianjie YU Rengui

（College of Information and Management Science， Henan Agricultural University，Zhengzhou Henan 450002）

Abstract： Newton interpolation method is one of the most important methods in numerical calculation. This paper introduced the theory of Newton's interpolation method， discussed its advantages and disadvantages and the relationship with Taylor's expansion.? Through a specific example， the realization and characteristics of Newton's interpolation method were shown intuitively by using MATLAB software.

Keywords： Newton interpolation method;Taylor expansion;Runge's phenomenon;MATLAB

1 研究背景

通常情況下，在科學(xué)研究和生產(chǎn)實踐中都用函數(shù)[y=fx]來表示某種內(nèi)在規(guī)律的數(shù)量關(guān)系。若[fx]在某個區(qū)間[a，b]上是存在的、連續(xù)的，但只能給出[a，b]上一系列點的函數(shù)值表（見表1）時，或者函數(shù)有解析表達式，但計算過于復(fù)雜、使用不方便只給出函數(shù)值表（如三角函數(shù)表、對數(shù)表等）時，為了研究函數(shù)的變化規(guī)律，往往需要求出不在表上點的函數(shù)值。因此，人們希望根據(jù)給定的函數(shù)表做一個既能反映函數(shù)[fx]的特性，又便于計算的簡單函數(shù)[Px]，用[Px]近似表示[fx]，即[Pxi=fxii=0，1，2，…，n]。這就引出了插值問題[1]。插值問題的幾何意義是：已知平面上n+1個不同的節(jié)點，要尋找一條次數(shù)不超過n的多項式曲線通過這些節(jié)點。

在我國古代，人們就開始運用插值法的思想來解決問題。隋代天文學(xué)家劉焯創(chuàng)立了用三次差內(nèi)插法來計算日月視差運動速度;唐代杰出天文學(xué)家張遂在《大衍歷》中提出了自變數(shù)不等間距的二次差內(nèi)插法用于編制天文數(shù)表。此外，17世紀(jì)之后，西方數(shù)學(xué)家拉格朗日和牛頓分別討論了非等距和等距的一般插值公式。目前，人們常用的插值方法有：拉格朗日插值、牛頓插值、埃爾米特插值、分段低次插值和三次樣條插值。

理論、試驗、計算是人類進行科學(xué)活動的三大方法，許多實際的科學(xué)與工程問題的解決都離不開科學(xué)計算。牛頓插值法是數(shù)值計算中較為重要的方法之一，在電力、機械、建筑科學(xué)與工程、環(huán)境科學(xué)、電信技術(shù)、計算機等學(xué)科中都有廣泛應(yīng)用。本文著重討論牛頓插值法理論以及MATLAB上機實驗。

2 牛頓插值法

在區(qū)間[a，b]上，函數(shù)[fx]關(guān)于一個節(jié)點[xi]的零階差商定義如式（1）所示，[fx]關(guān)于兩個節(jié)點[xi和xj]的一階差商定義見式（2）。一般地，k階差商就是k-1階差商的差商，稱式（3）為[fx]關(guān)于k+1個節(jié)點[x0，x1，…，xk]的k階差商[2]。具體可以按表2的格式有規(guī)律地計算差商。

[fxi=fxi]? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（1）

[fxi，xj=fxj-fxixj-xi]? ? ? ? ? ? ? ? ? ? ? ? ? ? （2）

[fx0，x1，…，xk=fx1，x2，…，xk-fx0，x1，…，xk-1xk-x0]? ? （3）

借助差商的定義，牛頓插值多項式可以表示為：

[Nnx=fx0ω0x+fx0，x1ω1x+fx0，x1，x2ω2x+…+fx0，x1，…，xnωnx]? ? ?（4）

牛頓插值多項式的余項公式可以表示為：

[Rnx=fx，x0，x1，…，xnωn+1x]? ? ? ? ? ? ? ?（5）

其中，[ω0x=1]，[ωkx=x-x0x-x1…x-xk-1][k=1，2，…，n+1]。對于[a，b]中的任一點[x]，則有[fx=Nnx+Rnx]。

假設(shè)由表1給出的插值問題，已經(jīng)得到牛頓插值多項式，如果所得插值多項式與原函數(shù)在某些點的誤差較大，達不到所要求的精度，一個很自然的想法就是希望通過增加節(jié)點以提高插值多項式的次數(shù)，從而提高精度。拉格朗日插值法的基函數(shù)簡單，且公式結(jié)構(gòu)緊湊，如果增加節(jié)點，所有基函數(shù)都需要重新進行計算，牛頓插值法就很好地克服了這一缺陷，在多項式中只需要增加一項即可。根據(jù)差商的性質(zhì)，對于在差值節(jié)點[x0，x1，…，xn]中任意位置增加或刪除一個插值節(jié)點的問題，可以設(shè)計牛頓插值多項式的承襲性算法。一般認(rèn)為，插值節(jié)點越多，插值多項式[Nnx]次數(shù)越高，逼近[fx]的效果就越好，但20世紀(jì)初龍格給出一個反例，隨著n的增大，高次多項式插值會出現(xiàn)震蕩現(xiàn)象，即龍格現(xiàn)象。

設(shè)[fx∈Cna，b]，一般稱[fx0]為[fx]在二重節(jié)點（[x0，x0]）處的一階差商，記作[fx0，x0=fx0]，稱[fx02！]為[fx]在三重節(jié)點[x0，x0，x0]處的二階差商，記作[fx0，x0，x0=fx02！]，依此類推，[fx]在n重節(jié)點[x0，x0，…，x0]處的n-1階差商記作[fx0，x0，…，x0=fx0n！]。由泰勒定理[3]可知，對[x∈a，b]，有

[fx=fx0+fx0x-x0+fx02！x-x02+…+fnx0n！x-x0n+Rnx]? ?（6）

其中，[Rnx=fn+1ξn+1！x-x0n+1]，[ξ]在[x]與[x0]之間。由[n]重節(jié)點差商的定義和泰勒定理可知：①通過計算點[x0]處的[n]階差商表，可以寫出該點處的[n]階泰勒展式;②通過計算點[x0]處的[n]階導(dǎo)數(shù)，可以寫出該點處的[n]階牛頓插值多項式;③牛頓插值余項與泰勒余項相等，即[fx]在[n]重節(jié)點[x0，x0，…，x0]處的牛頓插值多項式與[fx]在點[x0]處的泰勒展式相同。

3 牛頓插值法的MATLAB實現(xiàn)

例：設(shè)有函數(shù)[fx=11+x2]，在[-5，5]上取等距節(jié)點[xi=-5+ii=0，2，4，…，10]上寫出牛頓插值多項式，在[-5，5]上取等距節(jié)點[xi=-5+ii=0，2，4，…，10]上討論牛頓插值龍格現(xiàn)象。

解：在MATLAB命令窗口中輸入：

>>format rat

x0=-5：2：5; y0=1./（1+x0.^2）;

syms x; n=max（size（x0））;

y=y0（1）; disp（y）;

s=1; dx=y0;

for i=1：n-1

dx0=dx;

for j=1：n-i

dx（j）=（dx0（j+1）-dx0（j））/（x0（i+j）-x0（j））;

end

df=dx（1）; s=s*（x-x0（i））; y=y+s*df;

disp（y）;

end

N_n（x）=y

運行上述程序結(jié)果如下：

1/26

（2*x）/65 + 5/26

（2*x）/65 + （11*（x + 3）*（x + 5））/260 + 5/26

（2*x）/65 + （11*（x + 3）*（x + 5））/260 - （（x + 1）*（x + 3）*（x + 5））/65 + 5/26

（2*x）/65 + （11*（x + 3）*（x + 5））/260 - （（x + 1）*（x + 3）*（x + 5））/65 + （（x - 1）*（x + 1）*（x + 3）*（x + 5））/520 + 5/26

（2*x）/65 + （11*（x + 3）*（x + 5））/260 - （（x + 1）*（x + 3）*（x + 5））/65 + （（x - 1）*（x + 1）*（x + 3）*（x + 5））/520 + 5/26

N_n（x） =（2*x）/65 + （11*（x + 3）*（x + 5））/260 - （（x + 1）*（x + 3）*（x + 5））/65 + （（x - 1）*（x + 1）*（x + 3）*（x + 5））/520 + 5/26

在MATLAB命令窗口中輸入：

>>x0=-5：5; y0=1./（1+x0.^2）;

x1=-5：2：5; y1=1./（1+x1.^2）;

x2=-5：0.02：5; y2=1./（1+x2.^2）;

x=-5：0.05：5;

m=length（x）;

for k=1：m

tx=x（k）; n=length（x0）;

y=y0（1）; s=1; dx=y0;

for i=1：n-1

dx0=dx;

for j=1：n-i

dx（j）=（dx0（j+1）-dx0（j））/（x0（i+j）-x0（j））;

end

df=dx（1）;

s=s*（tx-x0（i））;

y=y+s*df;

end

yx0（k）=y;

end

for k=1：m

tx=x（k）; n=length（x1）;

y=y1（1）; s=1; dx=y1;

for i=1：n-1

dx1=dx;

for j=1：n-i

dx（j）=（dx1（j+1）-dx1（j））/（x1（i+j）-x1（j））;

end

df=dx（1）;

s=s*（tx-x1（i））;

y=y+s*df;

end

yx1（k）=y;

end

figure;

plot（x，yx0，'：g'，'LineWidth'，2）;

hold on; plot（x，yx1，'--k'，'LineWidth'，1.5）;

hold on; plot（x2，y2，'-b'，'LineWidth'，1）;

grid on; xlabel（'x軸'）; ylabel（'y軸'）;

legend（'n=10插值曲線'，'n=5插值曲線'，'原曲線'）;

運行上述程序結(jié)果如圖1所示。從圖中可以看出，相較于[N5x]在區(qū)間[-3，3]內(nèi)，[N10x]與[fx]擬合效果較好，在區(qū)間[-5，3]和[3，5]內(nèi)，[N10x]與[fx]偏離較遠，出現(xiàn)龍格現(xiàn)象。

4 結(jié)語

牛頓插值法可以通過計算差商表來求出，當(dāng)插值節(jié)點增加或減少時，不需要像拉格朗日插值法那樣重新計算。牛頓插值法具有承襲性，減少了計算量。當(dāng)插值節(jié)點增多時，牛頓插值會出現(xiàn)龍格現(xiàn)象。在給定重節(jié)點處的牛頓插值多項式與該點處的泰勒展式相同。

參考文獻：

[1]朱建新，李有法.數(shù)值計算方法[M].北京：高等教育出版社，2012.

[2]甄西豐.差商與牛頓插值多項式的承襲性算法[J].華中理工大學(xué)學(xué)報，2000（4）：35-38.

[3]韓云端.微積分概念解析[M].北京：高等教育出版社，2007.