周建杰 禹仁貴

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

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

中圖分類號(hào)：G642文獻(xiàn)標(biāo)識(shí)碼：A文章編號(hào)：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í)踐中都用函數(shù)[y=fx]來表示某種內(nèi)在規(guī)律的數(shù)量關(guān)系。若[fx]在某個(gè)區(qū)間[a，b]上是存在的、連續(xù)的，但只能給出[a，b]上一系列點(diǎn)的函數(shù)值表（見表1）時(shí)，或者函數(shù)有解析表達(dá)式，但計(jì)算過于復(fù)雜、使用不方便只給出函數(shù)值表（如三角函數(shù)表、對(duì)數(shù)表等）時(shí)，為了研究函數(shù)的變化規(guī)律，往往需要求出不在表上點(diǎn)的函數(shù)值。因此，人們希望根據(jù)給定的函數(shù)表做一個(gè)既能反映函數(shù)[fx]的特性，又便于計(jì)算的簡(jiǎn)單函數(shù)[Px]，用[Px]近似表示[fx]，即[Pxi=fxii=0，1，2，…，n]。這就引出了插值問題[1]。插值問題的幾何意義是：已知平面上n+1個(gè)不同的節(jié)點(diǎn)，要尋找一條次數(shù)不超過n的多項(xiàng)式曲線通過這些節(jié)點(diǎn)。

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

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

2 牛頓插值法

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

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

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

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

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

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

牛頓插值多項(xiàng)式的余項(xiàng)公式可以表示為：

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

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

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

設(shè)[fx∈Cna，b]，一般稱[fx0]為[fx]在二重節(jié)點(diǎn)（[x0，x0]）處的一階差商，記作[fx0，x0=fx0]，稱[fx02！]為[fx]在三重節(jié)點(diǎn)[x0，x0，x0]處的二階差商，記作[fx0，x0，x0=fx02！]，依此類推，[fx]在n重節(jié)點(diǎn)[x0，x0，…，x0]處的n-1階差商記作[fx0，x0，…，x0=fx0n！]。由泰勒定理[3]可知，對(duì)[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é)點(diǎn)差商的定義和泰勒定理可知：①通過計(jì)算點(diǎn)[x0]處的[n]階差商表，可以寫出該點(diǎn)處的[n]階泰勒展式;②通過計(jì)算點(diǎn)[x0]處的[n]階導(dǎo)數(shù)，可以寫出該點(diǎn)處的[n]階牛頓插值多項(xiàng)式;③牛頓插值余項(xiàng)與泰勒余項(xiàng)相等，即[fx]在[n]重節(jié)點(diǎn)[x0，x0，…，x0]處的牛頓插值多項(xiàng)式與[fx]在點(diǎn)[x0]處的泰勒展式相同。

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

例：設(shè)有函數(shù)[fx=11+x2]，在[-5，5]上取等距節(jié)點(diǎn)[xi=-5+ii=0，2，4，…，10]上寫出牛頓插值多項(xiàng)式，在[-5，5]上取等距節(jié)點(diǎn)[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

運(yùn)行上述程序結(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插值曲線'，'原曲線'）;

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

4 結(jié)語

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

參考文獻(xiàn)：

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

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

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