Metody numeryczne fizyki/Sposoby rozwiązywania układów równań różniczkowych zwyczajnych z pewnymi warunkami początkowymi

Sposoby rozwiązywania układów równań różniczkowych zwyczajnych z pewnymi warunkami początkowymi

---

Licencja

Autor: Mirosław Makowiecki Absolwent UMCS Fizyki Komputerowej Uniwersytetu Marii Curie-Skłodowskiej w Lublinie Email: miroslaw(kropka)makowiecki(małpa)gmail(kropka)pl Dotyczy: książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami. Użytkownika książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami nie zwalnia z odpowiedzialności prawnoautorskiej nieprzeczytanie warunków licencjonowania. Umowa prawna: Creative Commons: uznanie autorstwa, na tych samych warunkach, z możliwością obowiązywania dodatkowych ograniczeń. Autor tej książki dołożył wszelką staranność, aby informacje zawarte w książce były poprawne i najwyższej jakości, jednakże nie udzielana jest żadna gwarancja, czy też rękojma. Autor nie jest odpowiedzialny za wykorzystanie informacji zawarte w książce, nawet jeśli wywołaby jakąś szkodę, straty w zyskach, zastoju w prowadzeniu firmy, przedsiębiorstwa lub spółki bądź utraty informacji, niezależnie czy autor (a nawet Wikibooks) został powiadomiony o możliwości wystąpienie szkód. Informacje zawarte w książce mogą być wykorzystane tylko na własną odpowiedzialność.

---

---

Będziemy się zajmować układem m równań różniczkowych, których rozwiązaniem są funkcje y_i=y_i(x) dla i=1,...,m z warunkiem początkowym y_i(x₀)=y_i,0, przedstawionych jako:

${{dy_{i}} \over {dx}}=f_{i}(x,y_{1},...,y_{m})\;$

(7.1)

Jeśli będziemy dodatkowo przyjmować jeszcze oznaczenia, co do wektora y, f i y₀, których schemat jest:

$\mathbf {y} ={\begin{bmatrix}y_{1}\\\vdots \\y_{m}\end{bmatrix}}\;$

(7.2)

$\mathbf {f} ={\begin{bmatrix}f_{1}\\\vdots \\f_{m}\end{bmatrix}}\;$

(7.3)

$\mathbf {y} _{0}={\begin{bmatrix}y_{1,0}\\\vdots \\y_{m,0}\end{bmatrix}}\;$

(7.4)

Przy oznaczeniach (7.2), (7.3) i (7.4) układ równań różniczkowych (7.1) zapisujemy w postaci jednego równania różniczkowego:

${{d\mathbf {y} } \over {dx}}=\mathbf {f} (x,\mathbf {y} )\;$

(7.5)

Równanie (7.5) dla m=1 przechodzi w równanie skalarne. Przyjmijmy, że funkcja f(x,y) jest ciągła w zbiorze D={(x,y):x₀≤x≤b,y∈R^m}. Zakładamy, że funkcja f(x,y) wcześniej podanym zbiorze spełnia warunek Lipschitza względem zmiennej wektorowej Y, których dla dowolnych $\mathbf {y} ,{\tilde {\mathbf {y} }}\in \mathbf {R} ^{n}\;$, jeśli istnieje takie L dla x∈<x₀,b>, że jest spełniony warunek:

$||\mathbf {f} (x,\mathbf {y} )-\mathbf {f} (x,{\tilde {\mathbf {y} }})||\leq L||\mathbf {y} -{\tilde {\mathbf {y} }}||\;$

(7.6)

Można udowodnić, że przy powyższych założeniach istnieje przynajmniej jedno rozwiązanie zagadnienia (7.5). Rozwiązanie y_n będziemy wyznaczać, jeśli będziemy znali rozwiązania y_n.,y_n-1,...,y_n-k dla k≥0 przy zdefiniowanych punktach iteracji x_i=x₀+ih, w których będziemy wyznaczać funkcję f(x,y). Fakt takiego wyznaczania wyznaczania naszej funkcji wektorowej piszemy:

$\mathbf {A} _{n}(h,\mathbf {y} _{n-k},...,\mathbf {y} _{n};\mathbf {y} _{n+1})=0\;$

(7.7)

Ponieważ rozwiązania y_i wyznaczaliśmy z pewnym błędem, bo krok h jest skończony, a nie nieskończenie mały, to z prawej stronie w (7.7) powinna być pewna liczba ogólnie niezerowa przy zastąpieniu y_i przez ciąg wartości dokładnych Y_i(x). Nazwijmy ten wektor jako T_n, wtedy tą wielkość możemy przepisać do postaci liczby, która jest błędem wyznaczania wspomnianej wcześniej wielkości:

$\mathbf {T} _{n}=h^{p+1}\gamma _{p}+O(h^{p+2})\;$

(7.8)

We wzorze (7.8) wielkość γ_n jest wielkością stałą niezerową, a liczba "p" jest liczbą zwana rzędem dokładności lub też zwana też jest rzędem metody przybliżonej.

Rozwiązywanie układu równań różniczkowych metodą Eulera

Najbardziej prostą i łatwą metodą jest metoda Eulera, która jest rozwiązaniem równania różniczkowego wektorowego (7.5), a jego przedstawienie jest:

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+h\mathbf {f} (x_{n},\mathbf {y} _{n}){\mbox{ dla }}n\geq 1\;$

(7.9)

Będziemy się zajmować tutaj tą metoda tylko dla przypadku m=1. Weźmy sobie rozwiązanie dokładne rozwiązania (7.5), który wraz z błędem T_n przedstawiamy w formie:

$\mathbf {Y} _{n+1}=\mathbf {Y} _{n}+hf(x_{n},Y_{n})+T_{n}{\mbox{ gdzie }}T_{n}={{1} \over {2}}h^{2}Y^{''}(\xi _{n}){\mbox{ dla }}\xi _{n}\in (x_{n},x_{n+1})\;$

(7.10)

Przyjmując oszacowanie błędu ε_n=Y_n-y_n, wtedy odejmując równanie (7.9) od (7.10) otrzymujemy równanie, do którego wykorzystamy warunek Lipschitza, w postaci:

$\epsilon _{n+1}=\epsilon _{n}+h[f(x_{n},Y_{n})-f(x_{n},y_{n})]+T_{n}\leq (1+hL)|\epsilon _{n}|+|T_{n}|{\mbox{ gdzie }}|T_{n}|\leq {{1} \over {2}}h^{2}M_{2}=T\;$
${\mbox{ przy }}M_{2}=\max _{\langle x_{0},b\rangle }|Y^{''}(x)|\;$

(7.11)

Rozwińmy rekurencję (7.11) aż do samego parametru ε₀, który jest błędem wyznaczania wartości y₀, w rezultacie otrzymujemy błąd bezwzględny wartości elementu y_n:

$|\epsilon _{n}|\leq (1+hL)|\epsilon _{n-1}|+T\leq (1+hL)^{2}|\epsilon _{n-2}|+T+(1+hl)T\leq \cdot \;$
$\leq (1+hL)^{n}|\epsilon _{0}|+T+T(1+hL)^{2}+...+T(1+hL)^{n-1}=T{{(1+hL)^{n}-1} \over {hL}}+(1+hL)^{n}|\epsilon _{0}|\;$

(7.12)

Ponieważ zachodzi 1+u<e^u dla u>0, to równość (7.12) możemy napisać:

$|\epsilon _{n}|\leq {{T} \over {hL}}\left(e^{L(x_{n}-x_{0})}-1\right)+e^{L(x_{n}-x_{0})}|\epsilon _{0}|{\mbox{ dla }}nh=x_{n}-x_{0}\;$

(7.13)

Przyjmijmy oszacowanie ε₀=0 i wzór na T (7.11), wtedy wzór na błąd rozwiązania przybliżonego pisujemy w zależności od wskaźnika "n" i wartości początkowej x₀:

$|\epsilon _{n}|\leq {{1} \over {2}}{{M_{2}} \over {L}}h\left(e^{L(x_{n}-x_{0})}-1\right)\;$

(7.14)

Jeśli z oszacowania mamy h→0, to dla ustalonego punktu x (x-x₀=nh, n→∞, h→0) wartość y_n dąży do wartości dokładnej Y(x), a szybkość jego zbieżności do wartości dokładnej jest O(h).

Wprowadzenie do metod różnicowych

Niech "h" będzie krokiem całkowania, a a_i, b_i będą liczbami rzeczywistymi, wtedy możemy zbudować metodę różnicową, w którym za pomocą pewnych poprzednich wartości funkcji y_i i argumentów x_i możemy napisać następną wartość funkcji, przy pomocy równania:

$y_{n+1}=\sum _{i=1}^{k}a_{i}y_{n+1-i}+h\sum _{i=1}^{k}b_{i}f(x_{n+1-i},y_{n+1-i}){\mbox{ dla }}n\geq k-1\;$

(7.15)

Dla k≥1 metodę (7.15) możemy tą metodę zastosować do rozwiązywania wartości równania różniczkowego (7.1). Powyższą metodę stosujemy, gdy mamy elementy x₁,..,x_k-1, ale warunek początkowy gwarantuje nam tylko wartość funkcji w punkcie x₀, zatem aby wyznaczyć te nasze brakujące punkty należy posłużyć się innymi metodami. Gdy b₀=0, to powyższą metodę nazywamy wzorek ekstrapolacyjnym jawnym. Natomiast jeśli b₀≠0, to powyższy wzór przy pomocy wzoru (7.1) dla dokładnych obliczeń możemy przepisać w formie:

$0=Y_{n+1}-\sum _{i=1}^{k}a_{i}Y_{n+1-i}-h\sum _{i=1}^{k}b_{i}Y_{n+1-i}^{'}-T_{n}\;$

(7.16)

Wyrazy Y_n+1-i możemy rozłożyć w szereg Tayllora w otoczeniu punktu x_n+1-k, wtedy współczynniki A_i przy tych wyrazach przepisujemy w postaci:

$A_{0}=1-\sum _{i=1}^{k}a_{i}\;$

(7.17)

$A_{s}={{((n+1-0)-(n+1-k))^{s}} \over {s!}}-\sum _{i=1}^{k}{{((n+1-i)-(n+1-k))^{s}} \over {s!}}a_{i}-\;$
$-\sum _{i=1}^{k}b_{i}{{((n+1-i)-(n+1-k))^{s-1}} \over {(s-1)!}}=\;$
$={{k^{s}} \over {s!}}-\sum _{i=1}^{k}{{(k-i)^{s}} \over {s!}}a_{i}-\sum _{i=1}^{k}b_{i}{{(k-i)^{s-1}} \over {(s-1)!}}{\mbox{ dla }}s\geq 1{\mbox{ i }}0^{0}=1\;$

(7.18)

Jeśli A_i=0 i=0,1,...p, ale A_p+1≠0, to metoda (7.15) przy współczynnikach (7.17) i (7.18) nazywamy metodą różnicową rzędu p, ponadto błąd metody różnicowej wyznaczamy przy pomocy formuły:

$T_{n}=A_{p+1}h^{p+1}Y_{n-1-k}^{(p-1)}+O(h^{p+2})\;$

(7.19)

Metodę różnicową nazywamy zbieżną, gdy dla nh=x-x₀ mamy spełniony warunek:

$\lim _{\overset {h\rightarrow 0}{\underset {n\rightarrow \infty }{}}}y_{n}=Y(x)\;$

(7.20)

Twierdzenie
Warunkiem koniecznym zbieżności metody różnicowej jest, gdy współczynniki A₀ i A₁ są równe zero.

Analiza błędów

Błąd rozwiązania przybliżonego ${\tilde {y}}_{n}\;$ względem rozwiązania dokładnego Y_n określamy:

$\epsilon _{n}=Y_{n}-{\tilde {y}}\;$

(7.21)

Równanie iteracyjne określająca wartość dokładną iteracyjną przy błędzie iteracji T_n, określamy na całym przedziale, w której chcemy wyznaczyć rozwiązanie równania (7.1), jest napisana:

$Y_{n+1}=\sum _{i=1}^{k}a_{i}Y_{n+1-i}+h\sum _{i=1}^{k}b_{i}f(x_{n+1-i},Y_{n+1-i})+T_{n}\;$

(7.22)

Równanie iteracyjne na ciąg ${\tilde {y}}_{n}\;$ z błędem powstałym w wyniku przybliżeń na maszynie cyfrowej określamy:

${\tilde {y}}_{n+1}=\sum _{i=1}^{k}a_{i}{\tilde {y}}_{n+1-i}+h\sum _{i=1}^{k}b_{i}f(x_{n+1-i},{\tilde {y}}_{n+1-i})-\delta _{n}\;$

(7.23)

Wzory (7.22) i (7.23) odejmujemy je od siebie i wykorzystujemy wzór (7.21) i zakładając, że B_n=T_n+δ_n, a także przyjmując oznaczenie $\eta _{n+1-i}=f(x_{n+1-i},Y_{n+1-i})-f(x_{n+1-i},{\tilde {y}}_{n+1-i})\;$, w ten sposób otrzymujemy tożsamość:

$\epsilon _{n+1}=\sum _{i=1}^{k}a_{i}\epsilon _{n+1-i}+h\sum _{i=0}^{k}b_{i}\eta _{n+1-i}+B_{n}\;$

(7.24)

Będziemy wykorzystywać twierdzenie o wartości średniej, w ten sposób zmienną η_n+1-i możemy przepisać w postaci:

$\eta _{n+1-i}=f(x_{n+1-i},Y_{n+1-i})-f(x_{n+1-i},{\tilde {y}}_{n+1-i})=\;$
$=f_{y}(x_{n+1-i},{\tilde {y}}_{n+1-i},\epsilon _{n+1-i},\theta _{n+1-1}){\mbox{ dla }}\theta _{i}\in (0,1)\;$

(7.25)

Biorąc pod uwagę obliczenia (7.25), a także biorąc pod uwagę obliczenia (7.24) i przenosząc wyraz o wskaźniku i=0 w pierwszej sumie po prawej jego stronie na jej lewą stronę, w ten sposób otrzymujemy schemat dla n≥k-1:

$[I-hb_{0}f_{y}(x_{n+1},{\tilde {y}}_{n+1},\epsilon _{n+1},\theta _{n+1}]\epsilon _{n+1}=\;$
$=\sum _{i=1}^{k}[a_{i}I+hb_{i}f_{y}(x_{n+1-i},{\tilde {y}}_{n+1-i},\epsilon _{n+1-i},\theta _{n+1-i})]\epsilon _{n+1-i}+B_{n}\;$

(7.26)

Weźmy sobie, że pochodna funkcji f_y względem zmiennej y jest równa parametrowi Λ i jest równa pewnej stałej, tzn. Λ=f_y=const, która jest macierzą o stopniu m na m, bo mamy m równań (7.1) i m wartości x_i, dla i=1,..,m, co to można podstawić tak wytłumaczony wzór do (7.26) i otrzymać bardziej prostą do wyglądu równanie w porównaniu do poprzedniego dla n≥k-1:

$[\mathbf {I} -hb_{0}\mathbf {\Lambda } ]\mathbf {\epsilon } _{n+1}=\sum _{i=1}^{k}[a_{i}\mathbf {I} +hb_{i}\mathbf {\Lambda } ]\mathbf {\epsilon } _{n+1-i}+B\;$

(7.27)

Będziemy zakładać, że $\left(\sum _{i=0}^{k}b_{i}\right)det\Lambda \neq 0\;$, to rozwiązaniem równania (7.35) przy założeniu, że A₀ (7.17) jest równe zero, aby ciąg wartości (7.15) dążył do zera, przy n nieskończonym i malejącej różnicy h pomiędzy punktami x_i, czyli był zbieżny, co stąd wynika, że równanie (7.27) generuje ciąg ε_i:

$\epsilon _{n}=-{{1} \over {h\sum _{i=0}^{k}b_{i}}}\mathbf {\Lambda } ^{-1}\mathbf {B} \;$

(7.28)

Przyjmijmy natomiast λ=Λ, wtedy na podstawie tego równość (7.27) dla przypadku m=1 przyjmuje kształt:

$(1-h\lambda b_{0})\epsilon _{n+1}=\sum _{i=1}^{k}(a_{i}+h\lambda b_{i})\epsilon _{n+1-i}\;$

(7.29)

Rozwiązaniem równania (7.29) jest rozwiązaniem zależne od stałej μ i parametru n, i go przedstawiamy w formie:

$\epsilon _{n}=\mu ^{n}\;$

(7.30)

Rozwiązanie (7.30) podstawiamy do równania (7.29) i go obustronnie dzielimy przez μⁿ⁺¹ i mnożymy przez μ^k, w ten sposób otrzymujemy ostatecznie równanie:

$(1-h\lambda b_{0})\mu ^{k}=\sum _{i=1}^{k}(a_{i}+h\lambda b_{i})\mu ^{k-1}\;$

(7.31)

Równanie charakterystyczne (7.31) może mieć k różnych pierwiastków $\mu _{1}^{n},...,\mu _{k}^{n}\;$, wtedy rozwiązaniem równania (7.31), czyli ε_n jest rozwiązaniem w postaci:

$\epsilon _{n}^{1}=c_{1}\mu _{1}^{n}+...+c_{k}\mu _{k}^{n}\;$

(7.32)

a jeśli równanie (7.30) ma pierwiastki jednokrotne, a także różniące się od siebie, to rozwiązanie ε_n równania (7.30) możemy napisać:

$\epsilon _{n}=(c_{1}+c_{2}n+..+c_{s}n^{s-1})\mu _{0}^{n}+c_{s+1}\mu _{s+1}^{n}+...+c_{k}\mu _{k}^{n}\;$

(7.33)

Natomiast, jeśli (7.30) zawiera pierwiastki zespolone, to taki pierwiastek tego równania zapiszemy w postaci cosinusów i sinusów jako:

$\mu =|\mu |(\cos \theta \pm j\sin \theta )\Rightarrow \mu ^{n}=|\mu |^{n}(\cos n\theta \pm i\sin n\theta )\;$

(7.34)

to rozwiązaniem równania (7.29), gdy jego pierwiastki mogą być niektóre zespolone, piszemy przy pomocy obliczeń (7.34):

$\epsilon _{n}=(c_{1}\cos n\theta +c_{2}\sin n\theta )|\mu _{1}|^{n}+c_{1}\mu _{3}^{n}+...+\mu _{k}^{n}\;$

(7.35)

Zajmijmy się teraz rozwiązaniem równania (7.27), gdy m≠1, wtedy wzór na ε_n jest:

$\mathbf {\epsilon } _{n}=\mu ^{n}\mathbf {v} \;$

(7.36)

wtedy wzór (7.36) możemy podstawić do ostatnio wspomnianego równania, wtedy przy oznaczeniu:

$\mathbf {M} (\mu )=[\mathbf {I} -hb_{0}\mathbf {\Lambda } ]\mu ^{k}-\sum _{i=1}^{k}[a_{i}\mathbf {I} +hb_{i}\mathbf {\Lambda } ]\mu ^{k-i}\;$

(7.37)

wtedy przy powyższym oznaczeniu możemy powiedzieć:

$\mathbf {M} (\mu )\mathbf {v} =\mathbf {0} \;$

(7.38)

Aby rozwiązanie (7.38) miało niezerowe rozwiązania v, to musi zachodzi z wiadomości z algebry:

$det\mathbf {M} (\mu )=0\;$

(7.39)

Wzór (7.36) jest rozwiązaniem równania jednorodnego (7.38), przy definicji macierzy M(μ) (7.37), której wyznacznik się zeruje i on ma się w postaci funkcji zależnej od "km" i od μ_s, a także zależy od otrzymanej bazy wektorów naszego rozważanego równania jednorodnego, całkowite rozwiązanie jest jego kombinacją liniową:

$\mathbf {\epsilon } _{n}=\sum _{s=1}^{km}c_{s}\mu _{s}^{n}\mathbf {v} _{s}\;$

(7.40)

Dokonajmy przekształcenia macierzy (7.37), w taki sposób by do niego otrzymać macierz podobną, która jest macierzą trójkątną górną, w której na diagonali występują wartości własne, według ${\tilde {\mathbf {M} }}(\mu )=\mathbf {P} ^{-1}\mathbf {M} (\mu )\mathbf {P} \;$:

${\tilde {\mathbf {M} }}(\mu )=\mathbf {P} ^{-1}\mathbf {M} (\mu )\mathbf {P} =[(\mathbf {I} -hb_{0}{\tilde {\mathbf {\Lambda } }})\mu ^{k}-\sum _{i=1}^{n}(a_{i}\mathbf {I} +hb_{i}{\tilde {\mathbf {\Lambda } }})\mu ^{k-i}]\;$

(7.41)

Rozwiązaniem równania własnego (7.38) przy macierzy (7.41), gdy wyznacznik z macierzy (7.41) jest równy zero, wtedy równanie, z którego będziemy wyznaczać $\mu _{i}^{n}\;$, ma kształt dla j=1,...,m:

$(1-h\lambda _{j}b_{0})\mu ^{k}=\sum _{i=1}^{k}(a_{i}+h\lambda _{j}b_{i})\mu ^{k-1}\;$

(7.42)

Widzimy, że równanie (7.42) jest bardzo podobne do (7.30), dlatego go będziemy podobnie rozwiązywać.

Stabilność metody różnicowej i jej zbieżność do granicy

Zbadajmy rozwiązanie dla m=1, i napiszmy pochodną funkcji "y" względem wartości "y", którego to równanie przedstawimy wychodząc od (7.1) przy założeniu, że funkcja "f" jest wprost proporcjonalna do wartości funkcji "y":

$y^{'}=\lambda y\;$

(7.43)

wtedy równość (7.43) możemy podstawić do równania różnicowego (7.15), co w końcu do niego podstawimy y_k=μ^k otrzymując z niego wynikające równanie charakterystyczne w drugiej równości:

$(1-h\lambda b_{0})y_{n+1}=\sum _{i=1}^{k}(a_{i}+h\lambda b_{i})y_{n+1-i}\Rightarrow (1-h\lambda b_{0})\mu ^{k}=\sum _{i=1}^{k}(a_{i}+h\lambda b_{i})\mu ^{k-i}\;$

(7.44)

Równanie (7.44) jest bardzo podobne do równania (7.29), a więc tak samo będziemy rozwiązywać. Jego rozwiązanie jest:

$y_{n}=d_{1}\mu _{1}^{n}+...+d_{k}\mu _{k}^{n}\;$

(7.45)

Jeśli metoda różnicowa jest rzędu p, to pierwiastek charakterystyczny, który wynika z równania charakterystycznego wychodzącego od (7.44), możemy przedstawić w formie:

$\mu _{1}=e^{\lambda h}+O(h^{p+1})\Rightarrow \mu _{1}^{n}=(e^{\lambda h}+O(h^{p+1}))^{n}=e^{n\lambda h}+O(h^{p+1})\;$

(7.46)

Współczynniki d_i w (7.45) wynikają z "k" warunków początkowych. Wartość y_i dąży do y₀ przy h dążącym do zera dla i=1,...,k-1. Aby wyznaczyć współczynnik d₁(h) należy je obliczyć według:

$d_{1}=d_{1}(h)={{\begin{vmatrix}y_{0}&1&\cdot &1\\y_{1}&\mu _{2}&\cdot &\mu _{k}\\\cdots &\cdots &\cdots &\cdots \\y_{k-1}&\mu _{2}^{k-1}&\cdots &\mu _{k}^{k-1}\end{vmatrix}} \over {\begin{vmatrix}1&1&\cdot &1\\\mu _{1}&\mu _{2}&\cdot &\mu _{k}\\\cdots &\cdots &\cdots &\cdots \\\mu _{1}^{k-1}&\mu _{2}^{k-1}&\cdots &\mu _{k}^{k-1}\end{vmatrix}}}\;$

(7.47)

Gdy h→0, to współczynnik d_i powinien dąży do y₀, a pozostałe współczynnik przy tym samym h powinny dążyć do zera, co te dwa warunki możemy przepisać:

$\lim _{h\rightarrow 0}d_{i}(h)=y_{0}\;$

(7.48)

$\lim _{h\rightarrow 0}d_{i}(h)=0{\mbox{ dla }}j=2,3,..,k\;$

(7.49)

Zauważmy, że w rozwiązaniu (7.45) pierwszy wyraz dąży do rozwiązania dokładnego, patrząc na rozwiązanie (7.46), przedstawionych według Y(x)=y₀e^h(x-x₀) z warunkiem początkowym y(x₀)=y₀. Współczynniki μ₂,...,μ_k nazywamy pierwiastkami pasożytniczymi, a wyrazy w (7.45) oprócz pierwszego składnika nazywamy częścią pasożytniczą. Dla λ>0 i h w miarę małych możemy napisać μ₁=μ₁(h)>0, to wtedy błąd ε_n patrząc na (7.28) i na (7.32) przepisujemy w formie:

$\epsilon _{n}=c_{1}\mu _{1}^{n}+...+c_{k}\mu _{k}^{n}-{{B} \over {h\lambda \sum _{i=0}^{k}b_{i}}}\;$

(7.50)

Wielkość ma ściśle ustaloną i określoną wartość h>0, jeśli na pierwiastki równania charakterystycznego (7.44) nałożymy warunek |μ_i|≤1, to żaden że współczynników c_iμ_i nie będzie rósł, gdy "n" będzie dążyło do nieskończoności, a jeśli rozwiązanie dokładne będzie wzrastało, to nie należy oczekiwać |μ_i|≤1., wtedy $c_{1}\mu _{1}^{n}\;$ będzie wzrastało nieograniczenie. Jeśli położymy warunek |μ_i|≤|μ₁| gwarantuje by część pasożytnicza w (7.45) nie rosła szybciej niż pierwszy jego wyraz. Metodę różnicową nazywamy stabilną, gdy wartość bezwzględna każdego współczynnika μ_i w (7.45) była mniejsza lub równa jedynce, a gdy wartość bezwzględna μ₁ jest mniejsza lub równa wartości bezwzględej μ_i, to tą metodę nazywamy względnie stabilną. Powyższą metodę nazywamy absolutnie stabilną na przedziale (λβ), gdy spełniony jest warunek hλ∈(αβ). Metodę względnie stabilną na przedziale podobnie definiujemy. Przyjdźmy teraz dla przypadku m>1 i napiszmy równanie wynikającego z (7.1), które to przedstawimy w zależności od macierzy Λ wychodząc od równania (7.15):

$(\mathbf {I} -hb_{0}\mathbf {\Lambda } )\mathbf {y} _{n+1}=\sum _{i=1}^{k}(a_{i}\mathbf {I} +hb_{i}\mathbf {\Lambda } )\mathbf {y} _{n+1-i}\;$

(7.51)

Rozwiązaniem dokładnym równania (7.51) nazywamy wzór zależny od "h", który to przepiszemy w formie:

$\mathbf {Y} (x)=c_{1}e^{\lambda _{1}(x-x_{0})}\mathbf {v} _{1}+...+c_{m}e^{\lambda _{m}(x-x_{0})}\mathbf {v} _{m}\;$

(7.52)

Rozwiązaniem napisane dla poszczególnych punktów odległych od siebie o krok h przedstawiamy w zależności od współczynników μ_ij, która jest rozwiązaniem równania charakterystycznego wynikłego z (7.51):

$\mathbf {y} _{n}=\sum _{i=1}^{k}\sum _{j=1}^{m}d_{ij}\mu _{ij}^{n}\mathbf {v} _{ij}\;$

(7.53)

Błąd metody różnicowej (7.15) określamy na podstawie wzoru (7.27) przy wykorzystaniu (7.28) wzór na wyznaczanie dokładnego rozwiązania:

$\mathbf {\epsilon } _{n}=\sum _{i=1}^{k}\sum _{j=1}^{m}c_{ij}\mu _{ij}^{n}\mathbf {v} _{ij}-{{1} \over {h\sum _{i=0}^{k}b_{i}}}\mathbf {\Lambda } ^{-1}\mathbf {B} \;$

(7.54)

Metodę absolutnie stabilną i względnie stabilną określonych ε na współczynniki błędu określamy podobnie jak dla wyznaczania wartości y_i.

Twierdzenie
Aby metoda różnicowa byłaz bieżna względem metody różnicowej (7.15), to pierwiastki wielomianu:

$\rho (\mu )=\mu ^{k}-\sum _{i=1}^{k}a_{i}\mu ^{k-1}\;$

(7.55)

których wartość bezwzględna z μ powinna być mniejsza lub równa jeden, a jeśli ta wartość bezwzględna jest równa jeden, to wielomian (7.55) ma pierwiastki jednokrotne. Tą metodę spełniająca powyższe warunki nazywamy metodę stabilną w sensie Dahlquista.

Wyznaczmy teraz błąd rozwiązania dokładnego, która w przedziale <x₀,b> ma ciągłe pochodne, która jest rzędu p i rozwiązanie różnicowe jest zbieżne do rozwiązania dokładnego, wtedy jest prawdziwe oszacowanie:

$||{\tilde {\mathbf {y} }}_{n}-\mathbf {Y} _{n}||\leq \left\{{\tilde {c}}_{1}\theta +{\tilde {c}}_{2}(x_{n}-x_{0})\left({{\delta } \over {h}}+{\tilde {c}}_{2}||\mathbf {Y} ^{(p+1)}||h^{p}\right)\right\}e^{{\tilde {c}}_{4}(x_{n}-x_{0})}{\mbox{ dla }}n\geq k\;$

(7.56)

gdzie:

• $\theta =\max _{0\leq i\leq k-1}||\mathbf {y} _{i}-\mathbf {Y} _{i}||\;$,
• $\delta =\max _{n\geq k}||\delta _{n}||\;$,
• $||\mathbf {Y} ^{(p+1)}||=\max _{x_{0}\leq x\leq b}||\mathbf {Y} ^{(p+1)}(x)||\;$.

Wzory różnicowe Rungego-Kutta

Określmy metodę typu Rungego-Kutta, który jest pewnym wzorem analogicznym do wzoru różnicowego, określany jako:

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+\sum _{i=1}^{s}w_{i}\mathbf {K} _{i}\;$

(7.57)

a wzory na K_i określamy według dwóch wzorów:

$\mathbf {K} _{1}=h\mathbf {f} (x_{n},\mathbf {y} _{n})\;$

(7.58)

$\mathbf {K} _{i}=h\mathbf {f} \left(x_{n}+a_{i}h,\mathbf {y} _{n}+\sum _{j=1}^{i-1}b_{ij}\mathbf {K} _{j}\right){\mbox{ dla }}i>1\;$

(7.59)

• gdzie stałymi w w metodzie (7.57) i we wzorach (7.58) i (7.59) są w_ij, a_i, b_ij.

W przeciwieństwie do metody różnicowej w metodzie Rungego-Kutta wystarczy tylko warunek początkowy brzegowy. W prawej stronie metody (7.57) należy wykonywać liczenie funkcji f(x,y) s razy. Ale wykazano, że dla rzędu metody p=2,3,4 najmniejsza ich liczba obliczeń jest s=p, a dla p>4 ich liczba jest s>p. Rząd metody (7.57) wyznaczamy podobnie jak rząd metody różnicowej (7.15). Najbardziej praktyczne znaczenie ma metoda rzędu s=4, ale do wymaga wyczerpujących obliczeń. Bowiem nie wystarcza już porównanie współczynników przy "h", które te współczynniki zależą od f_i(x,y) i od ich pochodnych cząstkowych.

Metoda Ruungego-Kutta dla s=1 i w₁=1 jest metodą różnicową Eulera (7.9). Omówmy teraz bardzo prosty przypadek jak tylko możliwe, tzn. dla m=1 i s=2 Stosując twierdzenie Tayllora i rozkładając prawą i lewą stronę równania w naszej metodzie, otrzymujemy:

$Y_{n}+hY_{n}^{'}+{{h^{2}} \over {2}}Y_{n}^{''}+O(h^{3})=Y_{n}+h\left[w_{1}f+w_{2}\left(f+a_{2}h{{\partial f} \over {\partial x}}+b_{21}hf{{\partial f} \over {\partial y}}\right)\right]{\Bigg |}_{(x_{n},Y_{n})}+O(h^{3})\;$

(7.60)

gdzie:

• $Y_{n}^{'}=f|_{(x_{n},Y_{n})}\;$,  $Y_{n}^{''}=\left({{\partial f} \over {\partial x}}+f{{\partial f} \over {\partial y}}\right){\Bigg |}_{(x_{n},Y_{n})}\;$.

W powyższych obliczeniach korzystano ze wzoru (7.1). We wzorze możemy porównań współczynniki stojące przy hⁱ otrzymujemy stąd dwa wzory:

$w_{1}+w_{2}=1\;$

(7.61)

${{1} \over {2}}\left({{\partial f} \over {\partial x}}+f{{\partial f} \over {\partial y}}\right)=w_{2}a_{2}+w_{2}b_{21}f{{\partial f} \over {\partial y}}\;$

(7.62)

Chcemy by wzory w (7.6) niezależały od funkcji f(x,y), to należy przyjąć trzy bardzo ważne wzory:

$w_{1}+w_{2}=1\;$

(7.63)

$a_{2}w_{2}={{1} \over {2}}\;$

(7.64)

$b_{21}=a_{2}\;$

(7.65)

Jeśli będziemy przyjmować $Df={{\partial f} \over {\partial x}}+f{{\partial f} \over {\partial y}}\;$, $D^{2}f={{\partial ^{2}f} \over {\partial x^{2}}}+2f{{\partial ^{2}f} \over {\partial x\partial y}}+f^{2}{{\partial ^{2}f} \over {\partial y^{2}}}\;$, i będziemy przyjmować oszacowanie |f|≤M i $\left|{{\partial ^{i+j}f} \over {\partial x^{i}\partial y^{j}}}\right|\leq {{N} \over {M^{i-1}}}\;$ dla i+j≤2, to błąd metody Rundego-Kutty liczymy ze wzoru (7.8), a więc parametr tam występujący γ₂ przepisujemy w formie:

$\gamma _{2}=\left({{1} \over {6}}-{{1} \over {2}}a_{2}^{2}w_{2}\right)D^{2}f{\Bigg |}_{(x_{n},Y_{n})}+{{1} \over {6}}f_{y}Df{\Bigg |}_{(x_{n},Y_{n})}\leq \left[4\cdot \left({{1} \over {6}}-{{1} \over {2}}a_{2}^{2}w_{2}\right)+{{1} \over {3}}\right]MN\;$

(7.66)

Przy wyżej poczynionych uwagach metodę Rundego-Kutty przy poczynionych uwagach K₁=hf(x_n, y_n) i K₂=hf(x_n+a₂h,y_n+a₂K₁) przyjmujemy w postaci:

$y_{n+1}=y_{n}+\left(1-{{1} \over {2a_{2}}}\right)K_{1}+{{1} \over {2a_{2}}}K_{2}{\mbox{ dla }}a_{2}\neq 0\;$

(7.67)

Dla przypadku m>1 metoda Rundego-Kutty prowadzi do:

$\mathbf {Y} _{n}+h\mathbf {Y} _{n}^{''}+\mathbf {O} (h^{2})=\mathbf {Y} _{n}+h(\underbrace {w_{1}+w_{2}+...+w_{s}} _{1})\mathbf {f} (x_{n},\mathbf {Y} _{n})+\mathbf {O} (h^{2})\;$

(7.68)

Rozważania na temat stabilności metod Rungego-Kutty i jego rozwiązania

Weźmy sobie równanie (7.1), które tutaj przedstawimy w postaci y'=λy, zatem rozwiązanie równania (7.57) dla przypadku m=1 przedstawiamy:

$y_{n+1}=\left(1+h\lambda +O(h^{2})\right)y_{n}\;$

(7.69)

wtedy rozwiązaniem w metodzie Rundego-Kutty przedstawiamy w zależności od y₀, wychodząc od rozwiązania równania charakterystycznego, sposobem:

$\mu =1+h\lambda +O(h^{2})\Rightarrow y_{n}=\mu ^{n}=\left(1=h\lambda +O(h^{2})\right)^{n}y_{0}=\left(e^{\lambda h}+O(h^{2})\right]y_{0}=\;$
$=\left(e^{nh\lambda }+O(h^{2})\right)y_{0}=\left(e^{\lambda (x_{n}-x_{0})}+O(h^{2})\right)y_{0}\;$

(7.70)

Powyższe równanie dla h→0 dąży do równania dokładnego Y(x)=y₀e^λ(x_n-x₀). Już dla przypadku s=2, gdy rzędem metody jest p=s=2, to rozwiązaniem równania charakterystycznego nazywamy wzór poniżej, z którego możemy wyznaczyć ten sam wzór na rozwiązanie dokładne Y(x) jak uzyskano w postaci wynikania z punktu (7.70):

$\mu =1+h\lambda +{{1} \over {2}}h^{2}\lambda ^{2}\;$

(7.71)

Dla równania czwartego rzędu Rundego-Kutty układ jest określony przy pomocy parametrów w_i, a_ij, b_ij, który w sumie zawiera w sumie 13 niewiadomych i 11 równań, czyli nie są one jednoznacznie określone. W tabelce poniżej podano wzory dla metody czwartego rzędu. Przy metodzie Gilla obliczenia mają największą dokładność obliczeń. Wzór Ralstona daje największe osiągalne oszacowanie błędu γ_n. Rozwiązanie równania charakterystycznego dla metody Rundego-Kutty dla metody czwartego rzędu piszemy:

$\mu =1+h\lambda +{{1} \over {2}}h^{2}\lambda ^{2}+{{1} \over {6}}h^{3}\lambda ^{3}+{{1} \over {24}}h^{4}\lambda ^{4}\;$

(7.72)

	Wzór klasyczny	Wzór "trzech ósmych"	Wzór "Gilla"	Wzór Ralstona
w1	1/6	1/8	1/6	0,17476028
w2	1/3	3/8	$(1-1/{\sqrt {2}})/3\;$	-0,55148053
w3	1/3	3/8	$(1+1/{\sqrt {2}})/3\;$	1,20553547
w4	1/6	1/8	1/6	0.17118478
a2	1/2	1/3	1/2	0,4
a3	1/2	2/3	1/2	0,45573726
a4	1	1	1	1
b21	1/2	1/3	1/2	0,4
b31	0	-1/3	$({\sqrt {2}}-1)/2\;$	0,29697760
b32	1/2	1	$1-1/{\sqrt {2}}\;$	0,15875966
b41	0	1	0	0,21810038
b42	0	-1	$-1/{\sqrt {2}}\;$	-3,0505096470
b43	1	1	$1+1/{\sqrt {2}}\;$	3,83286432

Niektóre wzory różnicowe

Przedstawimy tutaj dziesięć wzorów różnicowych, w których wzory od 1-5 są wzorami jawnymi, a wzory 6-10 wzorami uwikłanymi, a oto te wzory:

Wzory różnicowe	p	$A_{p+1}\;$
$y_{n+1}=y_{n}+hy_{n}^{'}\;$	1	${{1} \over {2}}\;$
$y_{n+1}=y_{n}+{{h} \over {2}}(3y_{n}^{'}-y_{n-1}^{'})\;$	2	${{5} \over {12}}\;$
$y_{n+1}=y_{n}+{{h} \over {12}}(23y_{n}^{'}-16y_{n-1}^{'}+5y_{n-2}^{'})\;$	3	${{3} \over {8}}\;$
$y_{n+1}=y_{n}+{{h} \over {24}}(55y_{n}^{'}-59y_{n-1}^{'}+37y_{n-2}^{'}-9y_{n-3}^{'})\;$	4	${{251} \over {720}}\;$
$y_{n+1}=y_{n-3}+{{4} \over {3}}h(2y_{n}^{'}-y_{n-1}^{'}+2y_{n-2}^{'})\;$	4	${{14} \over {45}}\;$
$y_{n+1}=y_{n}+hy_{n+1}^{'}\;$	1	$-{{1} \over {2}}\;$
$y_{n+1}=y_{n}+(y_{n+1}^{'}+y_{n}^{'})\;$	2	$-{{1} \over {12}}\;$
$y_{n+1}=y_{n}+{{h} \over {12}}(5y_{n+1}^{'}+8y_{n}^{'}-y_{n-1}^{'})\;$	3	$-{{1} \over {24}}\;$
$y_{n+1}=y_{n}+{{h} \over {24}}(9y_{n+1}^{'}+19y_{n}^{'}-5y_{n-1}^{'}+y_{n-2}^{'})\;$	4	$-{{19} \over {720}}\;$
$y_{n+1}={{1} \over {8}}(9y_{n}-y_{n-2})+{{3} \over {8}}(y_{n+1}^{'}+2y_{n}^{'}-y_{n-1}^{'})\;$	4	$-{{1} \over {40}}\;$

Wzory 1-4 mają typ Adamsa-Bashfortha, a wzór 5 posiada typ Milne'a, a wzór 6-9 mają typ Adamsa-Moultona, a formuła 10 posiada typ Hamminga. Określmy teraz rząd metody zbieżnej, która wynosi k+2 dla k parzystego i dla k nieparzystego rząd jej wynosi k+1. Rozpatrujmy za Hammingiem wzór różnicowy, który przepisujemy do postaci:

$y_{n+1}=a_{1}y_{n}+a_{2}3y_{n-1}+a_{2}y_{n-2}+h(b_{0}y_{n+1}^{'}+b_{1}y_{n}^{'}+b_{2}y_{n-1}^{'})\;$

(7.73)

Jeśli zechcemy rozpatrzyć metodę czwartego rzędu równania różnicowego (7.73), co stąd otrzymamy sześć wzorów na współczynniki a₁, a₂, a₃, b₁, b₂, które zapiszemy od współczynnika a₂:

$a_{1}={{1} \over {8}}(9-a_{2})\;$

(7.74)

$a_{3}=-{{1} \over {8}}\left(1-a_{2}\right)\;$

(7.75)

$b_{0}={{1} \over {24}}(9-a_{2})\;$

(7.76)

$b_{1}={{1} \over {12}}(9-7a_{2})\;$

(7.77)

$b_{2}={{1} \over {24}}(-9+17a_{2})\;$

(7.78)

Błąd wyznaczania w metodzie różnicowej rozwiązania przybliżonego y_n względem rozwiązania dokładnego Y_n określamy metodą w zależności od a₂, wtedy:

$T_{n}=(-9+5a_{2})h^{5}y^{(5)}(\xi _{n})/360\;$

(7.79)

Wykorzystujemy podstawienie y_n=μⁿ, to wzór różnicowy (7.73) po wykorzystaniu tego ostatniego i wzory na pochodną y'_n=λ y_n, dochodzimy do wniosku że równanie charakterystyczne, które możemy przepisać w zależności od a₂:

$(1-h\lambda b_{0})\mu ^{3}=(a_{1}+h\lambda b_{1})\mu ^{2}+(a_{2}+h\lambda b_{2})\mu +a_{3}\;$

(7.80)

Dla $a_{2}\in \langle -0,6,1\rangle \;$ warunek stabilności jest spełniony przy h=0, ale już w przeciwnym przypadku musimy wybrać takie a₂ umieszczona w pobliżu środka przedziału, co wtedy wzór będzie stabilny dla tego przypadku, ale Hamming wziął przypadek a₂=0, wtedy wykorzystując równanie różnicowe (7.73) i wzory (7.74), (7.75), (7.76) i (7.77), to w takim przypadku równanie charakterystyczne (7.80), w którym przedziałem stabilności przy hλ piszemy (-2,6,0), rysujemy do postaci:

$\left({{3h\lambda } \over {8}}-1\right)\mu ^{3}=-\mu ^{2}\left({{9} \over {8}}+{{3h\lambda } \over {4}}\right)+\mu {{3h\lambda } \over {8}}+{{1} \over {8}}\;$

(7.81)

Weźmy sobie równanie różnicowe Newtona-Cotesa rozwiązujemy równanie róznicowe y'=f(x,y) przedziale zmienności <x_n-1,x_n+2>, wtedy rozważane równanie różnicowe przepisujemy do jego właściwej formy:

$y(x_{n+1})=y(x_{n-1})+\int _{x_{n-1}}^{x_{n}}f(x,y(x))dx\;$

(7.82)

Funkcja podcałkowa występująca w drugim składniku w równaniu różnicowym (7.82), w której występuje y(x), jest nie znana, więc zbudujmy odpowiedni wielomian L(x), który jest równoważny wielomianowi f(x,y(x)), który piszmy według wyglądu $_{L(x)=\sum _{i=0}^{l}f_{n-i}P_{i}(x)}\;$, wtedy równość (7.82) przepisujemy w formie:

$y_{n+1}=y_{n-1}+\sum _{i=0}^{l}f_{n-i}\int _{x_{n-1}}^{x_{n}}P_{i}(x)dx=y_{n-1}+h\sum _{i=0}b_{i}f_{n-i}\;$
${\mbox{ gdzie }}b_{i}=\int _{-j}^{l}P_{i}(x_{n}+th)dt{\mbox{ dla }}t={{x-x_{0}} \over {h}}\;$

(7.83)

Powyższy wzór nazywamy wzorem jawnym, ale już dla j=0 ten wzór przechodzi w wzór Adamsa-Bashfortha.

Zmienny etap całkowania

Zwykle dokonujemy całkowania przy stałym kroku całkowania, ale na maszynach cyfrowych stosuje zmienny krok całkowania. Na samym początku ustalamy taki krok całkowania, by wielkość błędu była mniejsza niż wielkość ε, to rozwiązanie równania całkowania powinno być żądaną dokładność w ściśle ustalonym przedziale całkowania, a jeśli błąd metody różnicowej jest większy od ustalonej wielkości, to ustalamy mniejszy krok całkowania. W metodzie różnicowej Rungego-Kutty stosuje tzw. zasadę Rungego, i ustalmy początkowo y_n=Y(x_n). Weźmy teraz przybliżone rozwiązanie metody różnicowej $y_{n+1}^{(1)}\;$, które stosujemy w punkcie x_n+1=x_n+h, dla obliczeń wykonanych dwukrotnie z krokiem h/2. Dla tego podwójnego kroku błąd metody różnicowej wykorzystując wzór (7.8) przepisujemy, gdzie p jest rzędem metody różnicowej, w formie:

$Y_{n+1}=y_{n+1}^{(1)}\simeq 2\gamma _{p}\left({{h} \over {2}}\right)^{p+1}\;$

(7.84)

A jeśli wykorzystamy krok całkowania h i będziemy obliczali wartość przybliżoną metody różnicowej $y_{n+1}^{(2)}\;$ ustalaną w punkcie x_n+1, wtedy błąd tejże metody jest napisany wzorem:

$Y_{n+1}-y_{n+1}^{(2)}\simeq \gamma _{p}h^{p+1}\;$

(7.85)

Jeśli będziemy odejmować wzory (7.84) i (7.85) od siebie, dalej i będziemy przyjmować fakt $_{Y_{n+1}-y_{n+1}^{(1)}\simeq \delta }\;$, wtedy na podstawie tego powiemy przy założeniu $_{\delta \simeq \gamma _{p}h^{p+1}{{1} \over {2^{p+1}}}}\;$:

$y_{n+1}^{(1)}-y_{n+1}^{(2)}=\gamma _{p}h^{p+1}\left(1-{{1} \over {2^{p}}}\right)=\delta \left(2^{p}-1\right)\Rightarrow \delta ={{y_{n+1}^{(1)}-y_{n}^{(2)}} \over {2^{p}-1}}\;$

(7.86)

Połączenie wzorów interpolacyjnych i ekstrapolacyjnych w metodach różnicowych

Weźmy sobie metodę różnicową (7.15). Jeśli w tej metodzie będziemy liczyli współczynniki A_p+1, to widzimy, że one dla b₀≠0 dla wzorów ektrapolacyjnych są większe niż dla wzorów interpolacyjnych. Dla tych formuł w metodzie ekstrapolacyjnych musimy na każdym etapie całkowania stosować wzory liniowe ze względu na y, podczas gdy we wzorach ekstrapolacji stosowane jest podstawienie do wzorów. Metodę iteracji prostej piszemy w zależności od ω_n w sposób:

$y_{n+1}^{(j+1)}=hb_{0}f(x_{n+1},y_{n+1}^{(j)})+\omega _{n}\;$

(7.87)

a w metodzie Newtona rysujemy:

$\mathbf {y} _{n+1}^{j+1}=\mathbf {y} _{n}^{(j)}+\left[\mathbf {I} -hb_{0}{{\partial f} \over {\partial \mathbf {y} }}(\mathbf {x} _{n},\mathbf {y} _{n}^{(j)})\right]\left(\mathbf {y} _{n}^{(j)}-hb_{0}f(\mathbf {x} _{n},\mathbf {y} _{n}^{(j)})-\omega _{n}\right)\;$

(7.88)

Wzory na współczynniki ω_n w (7.87) i (7.88) piszemy według schematu:

$\omega _{n}=\sum _{i=1}^{k}\left(a_{i}y_{n+1-i}+hb_{i}f(x_{n+1-i},y_{n+1-i})\right)\;$

(7.89)

Przy wykorzystaniu wzoru (7.89) wartość $\mathbf {y} _{n+1}^{(0)}\;$ należy dobrać w taki sposób by była ona bliska wartości $\mathbf {y} _{n+1}\;$, by można było wykonać najmniejszą liczbę iteracji, jak się okazuje, że najlepiej wartość $\mathbf {y} _{n+1}^{(0)}\;$ dobrać ze wzoru ekstrapolacyjnego, który nazwiemy wzorem wstępnym, a wzór iterpolacyjny, którego poprawiamy rozwiązanie $\mathbf {y} _{n+1}^{(0)}\;$ nazywamy wzorem korygującym. Jeśli połączymy wzory ekstrapolacyjne i interpolacyjne nazywamy metodą ekstrapolacyjno-interpolacyjną. Przy liczeniu kolejnych iteracji y_n stosuje się wzór wstępny i korygujący, a oto rozważmy przykład:

$y_{n+1}=y_{n}+hf(x_{n},y_{n})\;$

(7.90)

$y_{n+1}=y_{n}+{{h} \over {2}}\left[f(x_{n+1},y_{n+1})+f(x_{n},y_{n})\right]\;$

(7.91)

Z których wzór (7.90) nazywamy wzorem wstępnym, a drugi wzór (7.91) nazywamy wzorem korygującym, co w tych wzorach przedziały stabilności są (-2,0) i (-∞,0). Jeśli wykorzystamy wzór (7.87) i (7.88), w której to metodzie mamy jedną iterację prostą i korygującą, to zapisać będziemy mogli:

$y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})\;$

(7.92)

$y_{n+1}=y_{n+1}^{(1)}=y_{n}+{{h} \over {2}}\left[f(x_{n+1},y_{n+1}^{(0)})+f(x_{n},y_{n})\right]\;$

(7.93)

Metoda powyższa określona wzorami (7.92) i (7.93) jest metodą Rungego-Kutty drugiego rzędu, przedział stabilności tej metody jest (-2,0). A natomiast dokonamy jedną iterację wstepną i dwie iteracje korygujące, to otrzymamy wtedy:

$y_{n+1}^{(0)}=y_{n}+hf(x_{n},y_{n})\;$

(7.94)

$y_{n+1}^{(1)}=y_{n}+{{h} \over {2}}\left[f(x_{n+1},y_{n+1}^{(0)})+f(x_{n},y_{n})\right]\;$

(7.95)

$y_{n+1}=y_{n+1}^{(2)}=y_{n}+{{h} \over {2}}\left[f(x_{n+1},y_{n+1}^{(1)})+f(x_{n},y_{n})\right]\;$

(7.96)

Stabilność absolutna metody określonej trzema powyższymi metodami jest (-2,0).

Zastosowanie metody różnicowej Hamminga

Zastosujmy wzory 5 i 10 w ostatnio podanej tablicy, którego pierwszy jest wzorem wstępnym, a drugi wzorem korygującym, te formuły piszemy:

$y_{n+1}^{(0)}=y_{n-3}+{{4} \over {3}}\left(2y_{n}^{'}-y_{n-1}^{'}+2y_{n-2}^{'}\right)\;$

(7.97)

$y_{n+1}^{(j+1)}={{1} \over {8}}\left(9y_{n}-y_{n-2}\right)+{{3} \over {8}}h\left((y_{n+1}^{(j)})^{'}+2y_{n}^{'}-y_{n-1}^{'}\right)\;$

(7.98)

Równania przedstawiające równania na rozwiązania dokładne przy błędach $T_{n}^{0}\;$ dla wzoru wstępnego i przy $T_{n}\;$ dla wzory korygującego przedstawiamy:

$Y_{n+1}^{(0)}=Y_{n-3}+{{4} \over {3}}\left(2Y_{n}^{'}-Y_{n-1}^{'}+2Y_{n-2}^{'}\right)+T_{n}^{0}\;$

(7.99)

$Y_{n+1}^{(j+1)}={{1} \over {8}}\left(9Y_{n}-Y_{n-2}\right)+{{3} \over {8}}h\left((Y_{n+1}^{(j)})^{'}+2Y_{n}^{'}-Y_{n-1}^{'}\right)+T_{n}\;$

(7.100)

Błędy $T_{n}^{0}\;$, $T_{n}\;$ przedstawiamy jak nie trudno udowodnić wzorami:

$T_{n}^{0}={{14} \over {45}}h^{5}Y^{(5)}(\xi _{1}){\mbox{ dla }}\xi _{1}\in (x_{n-3},x_{n+1})\;$

(7.101)

$T_{n}=-{{1} \over {40}}h^{5}Y^{(5)}(\xi _{2}){\mbox{ dla }}\xi _{2}(x_{n-2},x_{n+1})\;$

(7.102)

Jeśli wprowadzimy oznaczenia ε_n=Y_n-y_n i ε'_n=Y'_n-y'_n, wtedy połączenie wzorów (7.97), (7.98) oraz (7.99), (7.100) daje nam w wyniku czego:

$y_{n+1}^{(0)}=Y_{n+1}-\epsilon _{n-3}-{{4} \over {3}}h(2\epsilon _{n}^{'}-\epsilon _{n-1}^{'}+2\epsilon _{n-2}^{'})-T_{n}^{0}\;$

(7.103)

$y_{n+1}=Y_{n+1}-{{1} \over {8}}(9\epsilon _{n}-\epsilon _{n-2})-{{3} \over {8}}(\epsilon _{n+1}^{'}+2\epsilon _{n}^{'}-\epsilon _{n-1}^{'})-T_{n}\;$

(7.104)

Jeśli połączymy dwa ostatnie wzory odejmując je od siebie w ostateczności dostajemy:

$y_{n+1}-y_{n+1}^{(0)}=-{{1} \over {8}}\left(9\epsilon _{n}-\epsilon _{n-2}\right)-h\left({{3} \over {8}}\epsilon _{n+1}^{'}+{{41} \over {12}}\epsilon _{n}^{'}-{{41} \over {24}}\epsilon _{n-1}^{'}+{{3} \over {8}}\epsilon _{n-2}^{'}\right)+T_{n}^{0}-T_{n}\;$

(7.105)

Jeśli założymy, że ε_i zmienia się wolno przy każdym kroku, a hε'_i jest małe od błędu metody i Y⁽⁵⁾ praktycznie się nie zmienia od ξ₁ i ξ₂, wtedy wykorzystując wzór (7.105) błędy dla rozwiązania przybliżonego względem dokładnego rozwiązania, tzn. wzory (7.101) i (7.102), piszemy:

$T_{n}^{0}={{112} \over {121}}(y_{n+1}-y_{n+1}^{(0)})\;$

(7.106)

$T_{n}=-{{9} \over {121}}(y_{n+1}-y_{n+1}^{(0)})\;$

(7.107)

Jeśli dodatkowo założymy, że różnica y_n+1 i $y_{n}^{(0)}\;$ jest mała z przejścia do przejścia, wtedy w (7.106) zamieniamy z n+1 na n. $T_{n}^{0}\;$ możemy zastąpić różnicą ${\overline {y}}_{n+1}^{(0)}-y_{n+1}^{(0)}\;$, a w tak powstałym równaniu $y_{n+1}^{(0)}\;$ wyznaczamy ze wzoru (7.97), a ${\overline {y}}_{n+1}^{(0)}\;$ piszemy ze wzoru z ostatnio powstałego wzoru po zastąpieniu $T_{n}^{0}\;$ przez różnice wcześniej powiedzianego, a ${\overline {y}}_{n+1}^{(0)}\;$ korygujemy wzorem korygującym przy pomocy wzoru (7.98) i wzory poniższego, którą nazywamy metodą Hamminga:

${\overline {y}}_{n+1}^{(0)}=y_{n+1}^{(0)}+{{112} \over {121}}(y_{n}-y_{n}^{(0)})\;$

(7.108)

Dla pierwszego przybliżenia $y_{n+1}^{(1)}\;$ bardzo dobrze przybliża rozwiązanie równania różnicowego y_n przy pierwszym zastosowaniu wzoru korygującego. To rozwiązanie możemy poprawić, bo potrafimy wyznaczyć błąd T_n. Na maszynach cyfrowych metodę Hamminga piszemy wedle schematów, w których pierwszy wzór jest wzorem wstępnym (7.97), a drugi modyfikujący (7.108), trzeci równanie różniczkowe, którego chcemy rozwiązać (7.1), w którym y zastępujemy przez ${\overline {y}}_{n+1}^{(0)}\;$, czwarty wzór korygujący zapisany według (7.98) dla j=0, w którym y_n+1 zastępowaliśmy przez ${\overline {y}}_{n+1}\;$, a piąty wzór końcowy jest to:

$y_{n+1}={\overline {y}}_{n+1}-{{9} \over {121}}({\overline {y}}_{n+1}-y_{n+1}^{(0)})\;$

(7.109)