LABORATORIUM

METODY NUMERYCZNE

Interpolacja metodą Lagrange'a i Newtona

Prowadzący: dr inż. Tomasz Kozłowski

Autor: Przemysław Wojciechowski

Treść zadania:

Dla funkcji f(x) = sin²(x) w punktach$\left\{ 0,\frac{\pi}{6},\frac{\pi}{2},\frac{3\pi}{4} \right\}$ wyprowadź ręcznie postacie wielomianu interpolacyjnego metodami wielomianów: Lagrange'a Oraz Newtona. Sprawdź, że
po uproszczeniach dają one identyczną postać wielomianu interpolacyjnego. Wyznacz następnie postacie obu wielomianów interpolacyjnych używając stosownych skryptów Matlaba. Do sprawozdania dołączyć należy zastosowany kod programu. Następnie sprawdź
i skomentuj we wnioskach błąd lokalny dla $x = \frac{\pi}{4}$. Wykreśl podane punkty {x_i,f(x₁)} przebieg funkcji f(x) oraz przebieg wyznaczonego wielomianu interpolacyjnego używając komendy plot. Wydruk wykresu (wraz z kodem) dołącz do sprawozdania.

Wstęp teoretyczny:

Wzór interpolacyjny Newtona ma następującą postać:

$$p_{k}\left( x \right) = \sum_{i = 0}^{n}{c_{i}\prod_{\begin{matrix} j = 0 \\ (j \neq i) \\ \end{matrix}}^{n}\left( x - x_{j} \right)}$$

Wzór ten można zapisać w wersji rozwiniętej:

p_n(x) = c₀ + c₁(x−x₀) + c₂(x−x₀)(x−x₁) + … + c_n(x−x₀)…(x−x_n − 1)

Współczynniki c_n obliczamy za pomocą ilorazów różnicowych np.:

$$c_{1} = f\lbrack x_{0},x_{1}\rbrack = \frac{f\left( x_{1} \right) - f\left( x_{0} \right)}{x_{1} - x_{0}}$$

Zauważmy, że inne uporządkowanie punktów (x₁,y₁) zupełnie zmienia składniki wzoru interpolacyjnego Newtona, choć nie narusza całej sumy. Wyraźmy wielomian L_n(x) jako sumę

$$L_{n}\left( x \right) = \sum_{i = 0}^{n}{y_{i}l_{i}\left( x \right)}$$

Po przekształceniach arytmetycznych otrzymujemy:

$$L_{n}\left( x \right) = \sum_{i = 0}^{n}{y_{i}\prod_{\begin{matrix} j = 0 \\ (j \neq i) \\ \end{matrix}}^{n}\frac{x - x_{j}}{x - x_{i}}}$$

Jeszcze inne wyrażenia wielomianu interpolacyjnego mają też swoje zalety i wady. Możemy np. szukać współczynników wielomianu przy potęgach zmiennej:

L_n(x) = λ₀(x)f(x₀) + λ₁(x)f(x₁) + … + λ_n(x)f(x_n)

Warunki interpolacyjne prowadzą do układu n + 1 równań liniowych wzglęgem tych współczynników:

$$\left\lbrack \begin{matrix} 1 & x_{0} & \ldots \\ \ldots & \ldots & \ldots \\ 1 & x_{n} & \ldots \\ \end{matrix}\text{\ \ \ \ }\begin{matrix} x_{0}^{n} \\ \ldots \\ x_{n}^{n} \\ \end{matrix} \right\rbrack \bullet \begin{bmatrix} a_{0} \\ \ldots \\ a_{n} \\ \end{bmatrix} = \begin{bmatrix} y_{0} \\ \ldots \\ y_{n} \\ \end{bmatrix}$$

Macierz tego układu nazywamy macierzą Vandermonde'a. Ta macierz bywa źle uwarunkowana i dlatego nie zaleca się stosowania tego układu. Zresztą koszt obliczeń byłby tu nadmiernie duży.

W obliczeniach numerycznych najbardziej użyteczny jest wzór interpolacyjny Newtona
w wersji zawierającej ilorazy różnicowe.

$$p\left( x \right) = \sum_{k = 0}^{n}{f\left\lbrack x_{1},x_{2},\ldots,x_{k} \right\rbrack\prod_{j = 0}^{k - 1}\left( x - x_{j} \right)}$$

Jego zaletą jest to, że dołączenie dodatkowych punktów (x₁,y₁) nie narusza obliczonych wcześniej współczynnikówc_j. Wartość tak wyrażonego wielomianu można łatwo obliczyć, stosując wariant schematu Hornera. Zaletą wzoru Lagrange'a jest natomiast niezależność wielomianuL_n(x) od rzędnych y_j, co przydaje się w rozważaniach analitycznych, np.
w konstrukcji kwadratur. Ten sam wzór może być wygodny wtedy, gdy dla ustalonych węzłów trzeba uwzględnić różne układy wielkości y_j, wynikających np. z pomiarów. Opracowano algorytm efektywnego obliczania wartości wielomianu za wzoru Lagrange'a (Werner [1984] i prace tam cytowane)

Przykład obliczeniowy:

Tab.1 Tabela z wartościami punktów

i	0	1	2	3
x	0	$$\frac{\mathbf{\pi}}{\mathbf{6}}$$	$$\frac{\mathbf{\pi}}{\mathbf{2}}$$	$$\frac{\mathbf{3}\mathbf{\pi}}{\mathbf{4}}$$
y	0	$$\frac{\mathbf{1}}{\mathbf{4}}$$	1	$$\frac{\mathbf{1}}{\mathbf{2}}$$

Obliczam wielomian metodą Lagrange'a:

L₃(x) = λ₀(x)f(x₀) + λ₁(x)f(x₁) + λ₃(x)f(x₃)

$$\lambda_{0}\left( x \right) = \frac{(x - x_{1})(x - x_{2})(x - x_{3})}{(x_{0} - x_{1})(x_{0} - x_{2})(x_{0} - x_{3})} = - 1,54x^{3} + 2,89x^{2} + 3,19x + 1,26$$

$$\lambda_{1}\left( x \right) = \frac{(x - x_{0})(x - x_{2})(x - x_{3})}{(x_{1} - x_{0})(x_{1} - x_{2})(x_{1} - x_{3})} = x^{3} - 3,93x^{2} + 3,7x$$

$$\lambda_{2}\left( x \right) = \frac{(x - x_{0})(x - x_{1})(x - x_{3})}{(x_{2} - x_{0})(x_{2} - x_{1})(x_{2} - x_{3})} = {- 0,775x}^{3} + 2,23x^{2} - 0,96x$$

$$\lambda_{3}\left( x \right) = \frac{(x - x_{0})(x - x_{1})(x - x_{2})}{(x_{3} - x_{0})(x_{3} - x_{1})(x_{3} - x_{2})} = {0,29x}^{3} - 0,6x^{2} + 0,24x$$

$$L_{3}\left( x \right) = \left( - 1,54x^{3} + 2,89x^{2} + 3,19x + 1,26 \right) \bullet 0 + \left( x^{3} - 3,93x^{2} + 3,7x \right) \bullet \frac{1}{4} =$$

$$+ \left( {- 0,775x}^{3} + 2,23x^{2} - 0,96x \right) \bullet 1 + \left( {0,29x}^{3} - 0,6x^{2} + 0,24x \right) \bullet \frac{1}{2} =$$

L3(x)=−0,38x^3+0,99x^2+0,085x

Obliczam wielomian metodą Newtona:

Q₃(x) = c₀ + c₁(x−x₀) + c₂(x−x₀)(x−x₁) + c₃(x−x₀)(x−x₁)(x−x₂)

c₀ = f(x₀) = 0

$$c_{1} = f\lbrack x_{0},x_{1}\rbrack = \frac{f\left( x_{1} \right) - f\left( x_{0} \right)}{x_{1} - x_{0}} = 0,48$$

$$c_{2} = f\left\lbrack x_{0},x_{1},x_{2} \right\rbrack = \frac{f\left\lbrack x_{1},x_{2} \right\rbrack - f\left\lbrack x_{0},x_{1} \right\rbrack}{x_{2} - x_{0}} = \frac{\frac{f\left( x_{2} \right) - f\left( x_{1} \right)}{x_{2} - x_{0}} - \frac{f\left( x_{1} \right) - f\left( x_{0} \right)}{x_{1} - x_{0}}}{x_{2} - x_{0}} = 0,15$$

$$c_{3} = f\left\lbrack x_{0},x_{1},x_{2},x_{3} \right\rbrack = \frac{f\left\lbrack x_{1},x_{2},x_{3} \right\rbrack - f\left\lbrack x_{0},x_{1},x_{2} \right\rbrack}{x_{3} - x_{0}} = \frac{\frac{f\left\lbrack x_{2},x_{3} \right\rbrack - f\left\lbrack x_{1},x_{2} \right\rbrack}{x_{3} - x_{1}} - 0,15}{x_{3} - x_{0}}$$

$$c_{3} = \frac{\frac{\frac{f\left( x_{3} \right) - f\left( x_{2} \right)}{x_{3} - x_{2}} - \frac{f\left( x_{2} \right) - f\left( x_{1} \right)}{x_{2} - x_{0}}}{x_{2} - x_{0}} - 0,15}{x_{3} - x_{0}} = - 0,365$$

Q₃(x) = c₀ + c₁(x−x₀) + c₂(x−x₀)(x−x₁) + c₃(x−x₀)(x−x₁)(x−x₂)

$$Q_{3}\left( x \right) = 0 + 0,48\left( x - 0 \right) + 0,15\left( x - 0 \right)\left( x - \frac{\pi}{6} \right) - 0,365\left( x - 0 \right)\left( x - \frac{\pi}{6} \right)\left( x - \frac{\pi}{2} \right)$$

Q3(x)=−0,37x^3+0,91x^2+0,01x

Kod źródłowy:

clc; clear all; close all; format short

%% Dane

xi = [0 pi/6 pi/2 3/4*pi];

q=inline('sin(x).*sin(x)');

yi =q(xi);

%% Interpolacja metoda Lagrange'a

n=length(xi);

ai=polyfit(xi,yi,n-1);

a=min(xi); b=max(xi);

x=linspace(a,b,100);

Ln=polyval(ai,x);

disp(sprintf('Wspolczynniki interpolacji Lagrangea'))

ai

%% Interpolacja metoda Newton

[cn,DD] = newtonp(xi,yi);

disp(sprintf('Wspolczynniki wielomianu Newtona'))

cn

Qn= polyval(cn,x);

%% Blad maksymalny

ksi=linspace(a,b);

f=inline('(-8*(cos(ksi)^2-sin(ksi)^2))/(4!)');

f=vectorize(f);

x=linspace(a,b);

fmax=max(abs(f(x)));

disp(sprintf('Blad maksymalny wynosi'))

fmax

%% Rysowanie wykresu

grid on

plot(xi,yi,'mo',x,q(x),'k',x,Ln,'.g',x,Qn,'b');

Skrypt obliczenia wielomianu Newtona:

function [cn,DD] = newtonp(x,y)

N = length(x) - 1;

DD = zeros(N+1,N+1);

DD(1:N+1,1)=y';

for k=2:N+1

for m=1:N+2-k

DD(m,k) = (DD(m+1,k-1) - DD(m,k-1))/(x(m+k-1)-x(m));

end

end

r = DD(1,:);

cn = r(N+1);

for k=N:-1:1

cn = [cn r(k)] - [0 cn*x(k)];

end

end

Wyniki obliczeń:

Tab.2 Wyniki obliczeń

Wielomian interpolacyjny wyliczony ręcznie metodą Lagrange'a	L3(x) = −0, 38x^3 + 0, 99x^2 + 0, 085x
Wielomian interpolacyjny wyliczony ręcznie metodą Newtona	Q3(x) = −0, 37x^3 + 0, 91x^2 + 0, 01x
Wielomian interpolacyjny wyliczony w programie Matlab	M3(x) = −0, 3778x^3 + 0, 9433x^2 + 0, 0872x

Wnioski:

Jak można zauważyć w tabeli 2 współczynnik przy najwyższej potędze dla wszystkich trzech wyliczonych wielomianów interpolacyjny jest bardzo zbliżony. Natomiast
w przypadku kolejnego współczynnika różnice są już zdecydowanie większe, lecz w oby przypadkach wielomianów obliczonych ręcznie są one prawie równo oddalone od wielomianu wyliczonego w programie Matlab. Ostatni współczynnik został najbliżej policzony metodą Lagrange'a w stosunku do wielomianu obliczonego w numerycznie. Różnice między wielomianami obliczonymi ręcznie a wielomianem obliczonym numerycznie wyniki z błędów zaokrągleń przy obliczeniach na kalkulatorze. Wartość przybliżenia przyjęta do celu uproszczenia tych obliczeń to trzy miejsca po przecinku.