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
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