Metoda Ekstrapolacji Richardsona

Przybliżanie, metoda najmniejszych kwadratów, wielomiany interpolacyjne i inne.
janusz47
Użytkownik
Użytkownik
Posty: 7910
Rejestracja: 18 mar 2009, o 16:24
Płeć: Mężczyzna
Podziękował: 30 razy
Pomógł: 1670 razy

Metoda Ekstrapolacji Richardsona

Post autor: janusz47 »

Metoda ekstrapolacji Richardsona

Rozwijamy funkcję odpowiedniej klasy w szereg Taylora w otoczeniu punktu \(\displaystyle{ x }\)

\(\displaystyle{ f(x+h) = f(x) + hf^{'}(x) + \frac{h^2}{2!}f^{''}(x) + \frac{h^3}{3!}f^{'''}(x)+\ \ ... \ \ (1)}\)

Obliczamy z równania \(\displaystyle{ (1) \ \ f^{'}(x) }\)

\(\displaystyle{ f^{'}(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f^{''}(x)- \frac{h^2}{6}f^{'''}(x)+ \ \ ... \ \ (2) }\)

Dla małych wartości \(\displaystyle{ h }\) równanie aproksymacyjne pierwszą pochodną \(\displaystyle{ (2) }\) możemy napisać w postaci

\(\displaystyle{ f^{'}(x) \approx \frac{f(x+h) -f(x)}{h} \ \ (3) }\)

Równanie \(\displaystyle{ (3) }\) nazywamy różnicą dzieloną w tył .

Błąd obcięcia aproksymacji \(\displaystyle{ f^{'}(x) }\) wynosi \(\displaystyle{ - \frac{h}{2}f^{''}(\zeta) }\) dla \(\displaystyle{ \zeta \in ( x-h), x). }\)

Jeśli w równaniu \(\displaystyle{ (1) }\) zastąpimy \(\displaystyle{ h }\) przez \(\displaystyle{ -h }\) i odejmiemy pierwsze równanie od drugiego to otrzymamy

najbardziej znaną różnicę aproksymującą pierwszą pochodną - tzw. różnicę centralną:

\(\displaystyle{ f^{'}(x) = \frac{f(x+h) -f(x-h)}{2h} - \frac{h^2}{3!}f^{'''}(x) - \frac{h^4}{5!}f^{(5)}- \ \ ... \ \ (3) }\)

Metoda ekstrapolacji Richardsona wykorzystuje równanie \(\displaystyle{ (3) }\)

Zapisujemy to równanie w postaci:

\(\displaystyle{ f^{'}(x) = ]\frac{f(x+h)-f(x-h)}{2h} +c_{2}h^2 + c_{4}h^{4} +... \ \ (4)}\)

gdzie współczynniki \(\displaystyle{ c_{2}, c_{4}, ... }\) zależą od \(\displaystyle{ f }\) i \(\displaystyle{ x.}\)

Dla ustalonej wartości \(\displaystyle{ x }\) definiujemy funkcję

\(\displaystyle{ \alpha(h) = \frac{f(x+h) -f(x-h)}{2h} \ \ (5) }\)

Z \(\displaystyle{ (4) }\)

\(\displaystyle{ f^{'}(x) = \alpha(h) + c_{2}h^2 + c_{4}h^{4} + \ \ ... \ \ (6) }\)

Pomysł Levvi Fry Richardsona na zmniejszenie błędu obcięcia polegał na wyeliminowaniu \(\displaystyle{ c_{2}h^2 }\) z równania \(\displaystyle{ (6).}\)

W tym celu zastępujemy w równaniu \(\displaystyle{ (6) \ \ h }\) przez \(\displaystyle{ \frac{h}{2}, }\)

\(\displaystyle{ f^{'}(x) = \alpha \left(\frac{h}{2}\right) + c_{2}\left(\frac{h}{2}\right)^2 + c_{4}\left(\frac{h}{2}\right)^4 + \ \ ... \ \ (7) }\)

i eliminujemy z równania \(\displaystyle{ (7) \ \ h^2 }\)

\(\displaystyle{ f^{'}(x) = \alpha \left (\frac{h}{2} \right) + \frac{1}{3}\left[ \alpha\left(\frac{h}{2}\right) -\alpha(h) \right] - \frac{1}{4}c_{4}h^4 + \ \ ...(8)}\)

Zaważmy, że dokładność aproksymacji pochodnej jest rzędu \(\displaystyle{ h^4. }\)

Możemy więc zapisać

\(\displaystyle{ f^{'}(x) \approx \alpha \left (\frac{h}{2} \right) + \frac{1}{3}\left[ \alpha\left(\frac{h}{2}\right) -\alpha(h) \right];}\)

Ta ekstrapolacyjna metoda może być kontynuowana do otrzymania dokładności rzędu \(\displaystyle{ h^6, \ \ h^8 }\) itd.

Przykład

Mamy daną funkcję \(\displaystyle{ f(x) = e^{x} }\) - aproksymujemy jej \(\displaystyle{ f^{'}(1,5) }\) z \(\displaystyle{ h = 0,1. }\)

\(\displaystyle{ \alpha(h) = \alpha(0,1) = \frac{f(1,6) - f(1,4)}{0,2} = \frac{e^{1,6} - e^{1,4}}{0,2} \approx 4, 488916, }\)

\(\displaystyle{ \alpha(h/2)= \alpha(0,05) = \frac{f(1,55) -f(1,45)}{01} = \frac{e^{1,55} - e^{1,45}}{0,1} \approx 4,483557.}\)

Na podstawie równania \(\displaystyle{ (8) }\)

\(\displaystyle{ f^{'}(1,5) \approx 4,483557 + \frac{1}{3}\left( 4, 483557 - 4,488916 \right )\approx 4,4817.}\)

W ogólnym przypadku obliczamy wartości

\(\displaystyle{ D_{n,1} = \alpha\left( \frac{h}{2^{n-1}}\right), \ \ n=1,2,\ \ ... }\) dla danego \(\displaystyle{ h>0.}\) tworząc kolumny z

wartościami \(\displaystyle{ D_{n,m}. }\)

Można udowodnić [ patrz np. ( *)], że zachodzi zależność rekurencyjna:

\(\displaystyle{ D_{n, m+1} \approx \frac{4^{m}}{4^{m} -1} D_{n,m} - \frac{1}{4^{m}-1} D_{n-1,m} \ \ (9) }\)

(*) David Kincaid Ward Cheney analiza numeryczna. WNT Warszawa 2006.


Program w MATLAB

Kod: Zaznacz cały

function eksrichardson(f,h,a,n)
% Aproksymacja pochodnej funkcji f dla x=a.
disp(' Tablica pochodnych ')
disp('____________________________________________________________________________________')
disp(' i      h          Di1          Di2          Di3             ...  ')
disp('_____________________________________________________________________________________')
D(1,1)=(feval(f,a+h)-feval(f,a-h))/(2*h);
fprintf('%2.0f %8.4f %12.4f\n',1,h,D(1,1));
for i=1:n-1
    h=h/2;
    D(i+1,1)=(feval(f,a+h)-feval(f,a-h))/(2*h);
    fprintf('%2.0f %8.4f %12.4f',i+1,h,D(i+1,1));
    for k=1:i
        D(i+1,k+1)=D(i+1,k)+(D(i+1,k)-D(i,k))/((4^k)-1);
        fprintf('%12.4f',D(i+1,k+1));
    end
    fprintf('\n');
end

[b]Przykład[/b]

function f=f1(x)
f=exp(x);

eksrichardson('f1',0.25,1,6)

Tablica pochodnych 
____________________________________________________________________________________________
 i      h          Di1          Di2          Di3             ...  
____________________________________________________________________________________________
 1   0.2500       2.7467
 2   0.1250       2.7254      2.7183
 3   0.0625       2.7201      2.7183      2.7183
 4   0.0313       2.7187      2.7183      2.7183      2.7183
 5   0.0156       2.7184      2.7183      2.7183      2.7183      2.7183
 6   0.0078       2.7183      2.7183      2.7183      2.7183      2.7183      2.7183

Przykład 2

Kod: Zaznacz cały

>> eksrichardson('f1',0.1,1.5,6)
 Tablica pochodnych 
____________________________________________________________________________________
 i      h          Di1          Di2          Di3             ...  
_____________________________________________________________________________________
 1   0.1000       4.4892
 2   0.0500       4.4836      4.4817
 3   0.0250       4.4822      4.4817      4.4817
 4   0.0125       4.4818      4.4817      4.4817      4.4817
 5   0.0063       4.4817      4.4817      4.4817      4.4817      4.4817
 6   0.0031       4.4817      4.4817      4.4817      4.4817      4.4817      4.4817
ODPOWIEDZ