Algorytm Neville’a

Algorytm Neville’a – algorytm zaproponowany przez angielskiego matematyka Erica Harolda Neville’a. Jest używany do wyznaczania wartości wielomianu interpolacyjnego (Lagrange’a i Newtona) w danym punkcie $x.$ Ideą jest wyznaczenie rozwiązania w krokach od pojedynczych węzłów do całego ich zbioru.

Biorąc pod uwagę zbiór danych punktów węzłowych $(x_i,y_i),i=1,2,\dots ,n,$ wielomian $P$ jest stopnia nie wyższego niż $n,$ a jego wartości w punktach węzłowych są takie same jak wartości interpolowanej funkcji:

$P(x_i)=f(x_i)=y_i,$     dla $i=1,\dots ,n+1.$

Definiujemy wielomiany interpolacyjne i ich wartości w ustalonym punkcie $x.$

• $p_i$ – wartość, w punkcie $x,$ wielomianu stopnia zerowego przechodzącego przez punkt $(x_i,y_i):$

$p_i=y_i=P_i(x)$     dla $i=0,\dots ,n.$

• $p_{i(i+1)}$ – wartość, w punkcie $x,$ wielomianu stopnia pierwszego przechodzącego przez punkty $(x_i,y_i)$ oraz $(x_{i+1},y_{i+1}):$

$p_{i(i+1)}=P_{i(i+1)}(x)$     dla $i=0,\dots ,n-1.$

• $p_{0\dots n}$ – wartość, w punkcie $x,$ wielomianu stopnia n-tego przechodzącego przez $(n+1)$ punktów $(x_i,y_i):$

$p_{0\dots n}=P_{0\dots n}(x)$     dla $i=0,\dots ,n.$

Wielomiany powyższego typu spełniają następującą własność rekurencyjną:

$p_{i(i+1)\dots (i+k)}=\frac{(x-x_{i+k})p_{i(i+1)\dots (i+k-1)}+(x_i-x)p_{(i+1)(i+2)\dots (i+k)}}{x_i-x_{i+k}},$

gdzie: $k$ odpowiada stopniowi wielomianu, $k=1,\dots ,n$ oraz $i=0,\dots ,n.$

Algorytm Neville’a polega na tym, że za pomocą powyższych wzorów konstruujemy tablicę symetryczną, która zawiera wartości wielomianu interpolacyjnego w ustalonym punkcie $x,$ dla $n=3:$

$x_0,y_0=p_0$
	$p_{01}$
$x_1,y_1=p_1$		$p_{012}$
	$p_{12}$		$p_{0123}$
$x_2,y_2=p_2$		$p_{123}$
	$p_{23}$
$x_3,y_3=p_3$

Kolejne elementy są obliczane rekurencyjnie na podstawie elementów poprzednich.

W praktyce algorytm Neville’a przedstawiamy w nieco innej wersji. Stosując oznaczenia:

$t_{(i+k),k}=p_{i(i+1),\dots ,(i+k)}$

Tablica przyjmuje postać:

$x_0,y_0=t_{00}$
	$t_{11}$
$x_1,y_1=t_{10}$		$t_{22}$
	$t_{21}$		$t_{33}$
$x_2,y_2=t_{20}$		$t_{32}$
	$t_{31}$
$x_3,y_3=t_{30}$

Ułatwia to komputerowe zaprogramowanie powyższej tablicy (jako tablicy dwuwymiarowej).

Otrzymujemy również wzór rekurencyjny w prostszej postaci:

$t_{ik}=\frac{(x-x_{i-k})t_{i,k-1}-(x-x_i)t_{i-1,k-1}}{x_i-x_{i-k}}$

gdzie: $1⩽k⩽i,$     dla $i=0,1,\dots ,n.$

Pseudokod:

       for i := 0 to n do
           t[i] = f[i]
       for i := i - 1 downto 0 do
           t[j]= t[j + 1] + (t[j + 1] - t[j]) * (x - x[i]) / (x[i] - x[j])

Szukaną wartość wielomianu interpolacyjnego otrzymujemy jako t[0].

Przykład

Dane mamy:

$i=0,1,2,3$

$x_0=0,$  $x_1=1,$  $x_2=3,$  $x_3=4$

$f(x_0)=1,$  $f(x_1)=3,$  $f(x_2)=2,$  $f(x_3)=1$

Wyliczymy wartość wielomianu interpolacyjnego w punkcie x=2. Korzystamy ze wzorów i konstruujemy tablicę:

$0,t_{00}=1$
	$t_{11}=5$
$1,t_{10}=3$		$t_{22}=10/3$
	$t_{21}=5/2$		$t_{33}=3$
$3,t_{20}=2$		$t_{32}=8/3$
	$t_{31}=3$
$4,t_{30}=1$

Zatem dla x=2 wartość wielomianu wynosi 3.

Uogólnienie algorytmu Neville’a

Mając zadane węzły $(x_i,f_i)$ do interpolacji można również użyć funkcji wymiernych:

$Z^{u,v}=\frac{P^{u,v}}{Q^{u,v}}=\frac{a_0+a_{1}x+\dots +a_u{x}^u}{b_0+b_{1}x+\dots +b_v{x}^v},$

spełniających warunki: $Z^{u,v}(x_i)=f_i,$     dla $i=0,1,\dots ,u+v.$

W przypadku interpolacji funkcjami wymiernymi można wyprowadzić uogólniony wzór Neville’a. Algorytm wygląda wówczas następująco:

$T_{i0}=f_i,T_{i,k-1}=0,T_{ik}=T_{i,k-1}+\frac{(T_{i,k-1}-T_{i-1,k-1})}{(\frac{x-x_{i-k}}{x-x_i}[1-\frac{T_{i,k-1}-T_{i-1,k-1}}{(T_{i,k-1}-T_{i-1,k-2}}]-1)}.$

Bibliografia

• J. Stoer, R. Bulirsch, Wstęp do analizy numerycznej, PWN, 1987.
• J.N. Lyness, C.B. Moler, Van Der Monde Systems and Numerical Differentiation, „Numerishe Mathematik” 8 (1966), s. 458–464