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

Z Wikibooks, biblioteki wolnych podręczników.
Przejdź do nawigacji Przejdź do wyszukiwania
Metody numeryczne fizyki
Metody numeryczne fizyki
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.makowiecki@gmail.com
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 oraz miejsca pochodzenia książki i jej jakikolwiek części, a także treści, teksty, tabele, wykresy, rysunki, wzory i inne elementy oraz ich części zawarte w książce, i tą książkę, nawet w postaci przerobionej nie można umieszczać w jakikolwiek formie na czasopismach naukowych, archiwach prac, itp.
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 yi=yi(x) dla i=1,...,m z warunkiem początkowym yi(x0)=yi,0, przedstawionych jako:

(7.1)

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

(7.2)
(7.3)
(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:

(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):x0≤x≤b,yRm}. Zakładamy, że funkcja f(x,y) wcześniej podanym zbiorze spełnia warunek Lipschitza względem zmiennej wektorowej Y, których dla dowolnych , jeśli istnieje takie L dla x∈<x0,b>, że jest spełniony warunek:

(7.6)

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

(7.7)

Ponieważ rozwiązania yi 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 yi przez ciąg wartości dokładnych Yi(x). Nazwijmy ten wektor jako Tn, wtedy tą wielkość możemy przepisać do postaci liczby, która jest błędem wyznaczania wspomnianej wcześniej wielkości:

(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[edytuj]

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

(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 Tn przedstawiamy w formie:

(7.10)

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


(7.11)

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


(7.12)

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

(7.13)

Przyjmijmy oszacowanie ε0=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 x0:

(7.14)

Jeśli z oszacowania mamy h→0, to dla ustalonego punktu x (x-x0=nh, n→∞, h→0) wartość yn 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[edytuj]

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

(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 x1,..,xk-1, ale warunek początkowy gwarantuje nam tylko wartość funkcji w punkcie x0, zatem aby wyznaczyć te nasze brakujące punkty należy posłużyć się innymi metodami. Gdy b0=0, to powyższą metodę nazywamy wzorek ekstrapolacyjnym jawnym. Natomiast jeśli b0≠0, to powyższy wzór przy pomocy wzoru (7.1) dla dokładnych obliczeń możemy przepisać w formie:

(7.16)

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

(7.17)


(7.18)

Jeśli Ai=0 i=0,1,...p, ale Ap+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:

(7.19)

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

(7.20)
Twierdzenie
Warunkiem koniecznym zbieżności metody różnicowej jest, gdy współczynniki A0 i A1 są równe zero.

Analiza błędów[edytuj]

Błąd rozwiązania przybliżonego względem rozwiązania dokładnego Yn określamy:

(7.21)

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

(7.22)

Równanie iteracyjne na ciąg z błędem powstałym w wyniku przybliżeń na maszynie cyfrowej określamy:

(7.23)

Wzory (7.22) i (7.23) odejmujemy je od siebie i wykorzystujemy wzór (7.21) i zakładając, że Bn=Tnn, a także przyjmując oznaczenie , w ten sposób otrzymujemy tożsamość:

(7.24)

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


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


(7.26)

Weźmy sobie, że pochodna funkcji fy względem zmiennej y jest równa parametrowi Λ i jest równa pewnej stałej, tzn. Λ=fy=const, która jest macierzą o stopniu m na m, bo mamy m równań (7.1) i m wartości xi, 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:

(7.27)

Będziemy zakładać, że , to rozwiązaniem równania (7.35) przy założeniu, że A0 (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 xi, czyli był zbieżny, co stąd wynika, że równanie (7.27) generuje ciąg εi:

(7.28)

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

(7.29)

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

(7.30)

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

(7.31)

Równanie charakterystyczne (7.31) może mieć k różnych pierwiastków , wtedy rozwiązaniem równania (7.31), czyli εn jest rozwiązaniem w postaci:

(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ć:

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

(7.34)

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

(7.35)

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

(7.36)

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

(7.37)

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

(7.38)

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

(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ą:

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

(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ć , ma kształt dla j=1,...,m:

(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[edytuj]

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

(7.43)

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

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

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

(7.46)

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

(7.47)

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

(7.48)
(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)=y0eh(x-x0) z warunkiem początkowym y(x0)=y0. Współczynniki μ2,...,μ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ć μ11(h)>0, to wtedy błąd εn patrząc na (7.28) i na (7.32) przepisujemy w formie:

(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 ciμ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 będzie wzrastało nieograniczenie. Jeśli położymy warunek |μi|≤|μ1| 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 μ1 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):

(7.51)

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

(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):

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

(7.54)

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

Twierdzenie
Aby metoda różnicowa byłaz bieżna względem metody różnicowej (7.15), to pierwiastki wielomianu:
(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 <x0,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:

(7.56)

gdzie:

  • ,
  • ,
  • .

Wzory różnicowe Rungego-Kutta[edytuj]

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

(7.57)

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

(7.58)
(7.59)

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 fi(x,y) i od ich pochodnych cząstkowych.

Metoda Ruungego-Kutta dla s=1 i w1=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:

(7.60)

gdzie:

  • .

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

(7.61)
(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:

(7.63)
(7.64)
(7.65)

Jeśli będziemy przyjmować, i będziemy przyjmować oszacowanie |f|≤M i dla i+j≤2, to błąd metody Rundego-Kutty liczymy ze wzoru (7.8), a więc parametr tam występujący γ2 przepisujemy w formie:

(7.66)

Przy wyżej poczynionych uwagach metodę Rundego-Kutty przy poczynionych uwagach K1=hf(xn, yn) i K2=hf(xn+a2h,yn+a2K1) przyjmujemy w postaci:

(7.67)

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

(7.68)

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

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

(7.69)

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


(7.70)

Powyższe równanie dla h→0 dąży do równania dokładnego Y(x)=y0eλ(xn-x0). 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):

(7.71)

Dla równania czwartego rzędu Rundego-Kutty układ jest określony przy pomocy parametrów wi, aij, bij, 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:

(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
-0,55148053
w3
1/3
3/8
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
0,29697760
b32
1/2
1
0,15875966
b41
0
1
0
0,21810038
b42
0
-1
-3,0505096470
b43
1
1
3,83286432

Niektóre wzory różnicowe[edytuj]

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
1
2
3
4
4
1
2
3
4
4

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:

(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 a1, a2, a3, b1, b2, które zapiszemy od współczynnika a2:

(7.74)
(7.75)
(7.76)
(7.77)
(7.78)

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

(7.79)

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

(7.80)

Dla warunek stabilności jest spełniony przy h=0, ale już w przeciwnym przypadku musimy wybrać takie a2 umieszczona w pobliżu środka przedziału, co wtedy wzór będzie stabilny dla tego przypadku, ale Hamming wziął przypadek a2=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:

(7.81)

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

(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 , wtedy równość (7.82) przepisujemy w formie:


(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[edytuj]

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 yn=Y(xn). Weźmy teraz przybliżone rozwiązanie metody różnicowej , które stosujemy w punkcie xn+1=xn+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:

(7.84)

A jeśli wykorzystamy krok całkowania h i będziemy obliczali wartość przybliżoną metody różnicowej ustalaną w punkcie xn+1, wtedy błąd tejże metody jest napisany wzorem:

(7.85)

Jeśli będziemy odejmować wzory (7.84) i (7.85) od siebie, dalej i będziemy przyjmować fakt , wtedy na podstawie tego powiemy przy założeniu :

(7.86)

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

Weźmy sobie metodę różnicową (7.15). Jeśli w tej metodzie będziemy liczyli współczynniki Ap+1, to widzimy, że one dla b0≠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:

(7.87)

a w metodzie Newtona rysujemy:

(7.88)

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

(7.89)

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

(7.90)
(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:

(7.92)
(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:

(7.94)
(7.95)
(7.96)

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

Zastosowanie metody różnicowej Hamminga[edytuj]

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:

(7.97)
(7.98)

Równania przedstawiające równania na rozwiązania dokładne przy błędach dla wzoru wstępnego i przy dla wzory korygującego przedstawiamy:

(7.99)
(7.100)

Błędy , przedstawiamy jak nie trudno udowodnić wzorami:

(7.101)
(7.102)

Jeśli wprowadzimy oznaczenia εn=Yn-yn 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:

(7.103)
(7.104)

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

(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(5) praktycznie się nie zmienia od ξ1 i ξ2, 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:

(7.106)
(7.107)

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

(7.108)

Dla pierwszego przybliżenia bardzo dobrze przybliża rozwiązanie równania różnicowego yn przy pierwszym zastosowaniu wzoru korygującego. To rozwiązanie możemy poprawić, bo potrafimy wyznaczyć błąd Tn. 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 , czwarty wzór korygujący zapisany według (7.98) dla j=0, w którym yn+1 zastępowaliśmy przez , a piąty wzór końcowy jest to:

(7.109)