Wiederholung: Landau-Notation und Taylor-Entwicklung

Landau-Notation:  Seien f,g:]0,+[Rn Abbildungen.
Man schreibt f=O(g), falls c>0δ>0x]0,δ[f(x)cg(x).
Man schreibt f=o(g), falls ε>0δ>0x]0,δ[f(x)εg(x).

Beispiel: f=O(1) gilt genau dann, wenn f in einer δ-Umgebung von 0 beschränkt ist.
f=o(1) ist äquivalent zu limx0f(x)=0.
f=o(x) ist äquivalent zu f~=o(1) mit f(x)=xf~(x).

Satz (Taylor-Entwicklung):
Seien UR ein Intervall und f:URRn in x0U (m+1)-fach stetig differenzierbar.
Dann gilt f(x0+h)=f(x0)+k=1m1k!f(k)(x0)hk+rm(x0,h) mit
rm(x0,h)=1(m+1)!f(m+1)(y)hm+1 für ein yx0,x0+h¯, d. h.
rm(x0,h)=O(hm+1). Es gilt auch rm(x0,h)=o(hm).

Motivation, Beispiele

Bemerkung: Gegeben seien eine Funktion f:R×RR, t0R und u0R. Gesucht ist eine differenzierbare Funktion u=u(t):RR, sodass u(t)=f(t,u(t)) für tt0. Dieses Problem heißt Anfangswertproblem.

Beispiel: Sei t die Zeit und P(t) eine Population. Für die Zunahme ΔP:=P(t+Δt)P(t) in der Zeit Δt soll ΔPαP(t)Δt mit α>0 gelten. Für Δt0 erhält man die DGL dPdt=αP(t). Sie hat die allgemeine Lösung P(t)=ceαt mit c beliebig (exponentielles Wachstum). Ist ein Anfangswert P0=P(t0) gegeben, so bestimmt sich c durch c=P0eαt0, d. h. die partikuläre Lösung ist P(t)=P0eα(tt0).

Beispiel: Die DGL dPdt=λP(KP) mit λ,K>0 modelliert logistisches Wachstum. Zum Beispiel gilt für PK, dass dPdt=0, d. h. P ändert sich nicht. Die DGL hat die Lösung P(t)=K1+KP01eλKt.

Beispiel: Eine DGL, mit der das aktuelle Bevölkerungswachstum beschrieben werden kann, lautet dPdt=αP(t)β mit α>0, β>1.

Beispiel: Wird die Menge einer radioaktiven Substanz durch u=u(t) beschrieben, so modelliert die DGL du=λudt, λ>0 den Zerfall der Substanz aufgrund der Radioaktivität. Für t0=0 lautet eine Lösung u(t)=u0eλt. Die Halbwertszeit ist die Zeit, in der sich die Menge der Substanz halbiert. Sie ist unabhängig von der aktuellen Menge und beträgt τ=ln(2)λ.

Theoretische Grundlagen

Anfangswertproblem:  Seien URn offen (Zustandsraum), fC(R×U,Rn), u0U, IR und t0I. Gesucht ist eine Funktion u=(u1,,un)tC1(I,U) mit u(t)=f(t,u(t)) für tI und u(t0)=u0. Dieses Problem heißt Anfangswertproblem (AWP).

Beispiel: Im Räuber-Beute-Modell wird mit y1(t) bzw. y2(t) die Population der Beute- bzw. Raubtiere bezeichnet. Die DGLs y1(t)=αy1(t)(1y2(t)) und y2(t)=βy2(t)(y1(t)1) modellieren dann den zeitlichen Verlauf der Populationen.

Existenz und Eindeutigkeit der Lösung des Anfangswertproblems

Satz (Peano): Seien f auf einem kompakten Rechteck
R:={(t,u)R×U||tt0|a,uu0b} stetig,
μ:=max(t,u)Rf(t,u)< und α:=min(a,bμ).
Dann hat das Anfangswertproblem auf [t0α,t0+α] mindestens eine Lösung.

Satz (Picard-Lindelöf): Sei zusätzlich f in R im zweiten Argument Lipschitz-stetig, d. h. f(t,w)f(t,w~)Lww~ für alle (t,w),(t,w~)R.
Dann existiert für U=Rn genau eine Lösung uC1([t0α,t0+α],Rn).

Satz (Banachscher Fixpunktsatz): Seien (X,) ein Banachraum und DX eine abgeschlossene Teilmenge mit D. Sei außerdem T:DX eine Abbildung mit T(D)D und 0<c<1v,v~DTvTv~cvv~. Dann gibt es genau ein uD, sodass Tu=u.

Behandlung von Anfangswertproblemen höherer Ordnung

Bemerkung: Ein Anfangswertproblem höherer Ordnung ist ein Anfangswertproblem der Form y(m)(t)=f(t,y(t),y(t),,y(m1)(t)) mit y(i)(t0)=y0,i für i=0,,m1.
Es kann in ein System 1. Ordnung umgeformt werden, indem z1(t)=y(t), z2(t)=y(t), …, zm(t)=y(m1)(t) gesetzt wird. Damit ist z=(z1,z2,,zm1,zm)t
=(z2,z3,,zm,f(t,z1,z2,,zm))t ein System 1. Ordnung mit der Anfangsbedingung
z(t0)=(y(t0),y(t0),,y(m1)(t0))t=(y0,0,y0,1,,y0,m1)t.

Beispiel: Die elastische Schwingung eines fest eingespannten Federpendels, an dem ein Körper mit Masse m hängt, kann durch die DGL my(t)+ry(t)+D(y(t))=g(t) beschrieben werden, wenn y(t) die Auslenkung darstellt und y(0) und y(0) gegeben sind. Umgeformt nach y ergibt dies y=1m(gD(y)ry). Mit z1=y und z2=y ist z=(z1,z2)t=(z2,1m(gD(z1)rz2))t ein System 1. Ordnung mit Anfangsbedingung z(0)=(y0,0,y0,1)t.

Lösung durch Trennung der Variablen

Bemerkung: Eine DGL hat trennbare Veränderliche, falls sie die Form y(t)=f(t)g(y) mit y(t0)=y0 besitzt. In diesem Fall kann sie mit der Gleichung 1g(y)dy=f(t)dt und anschließendem Integrieren, also y0y1g(z)dz=t0tf(s)ds, gelöst werden, indem nach y(t) umgeformt und die Integrationskonstante mit der Anfangsbedingung berechnet wird.

Satz (Korrektheit der Trennung der Veränderlichen): Seien fC(It,R), gC(Iy,R) und t0 bzw. y0 seien aus dem Inneren von It bzw. Iy. In diesem Fall ist die obige DGL mit dem eben beschriebenen Algorithmus in einer Umgebung von t0 eindeutig lösbar.

Spezielle Typen von DGL 1. Ordnung

autonom:  Eine DGL u(t)=f(t,u(t)) heißt autonom, falls u(t)=f(u(t)).

linear:  Eine DGL u(t)=f(t,u(t)) heißt linear, falls u(t)=A(t)u(t)+b(t) mit
AC(I,Rn×n) und bC(I,Rn).
Eine lineare DGL heißt homogen, falls b0, sonst heißt sie inhomogen/affin.

Satz (eindeutige Lösbarkeit linearer DGLs):
Sei u(t)=A(t)u(t)+b(t) eine lineare DGL mit AC(I,Rn×n)L(R,Rn×n).
Dann hat das Anfangswertproblem genau eine Lösung in C1(I,Rn).

Satz (Lösungen linearer DGLs): Unter den Voraussetzungen von eben gilt:

  • Die Lösungen der homogenen DGL u(t)=A(t)u(t) bilden einen n-dimensionalen Unterraum VC1(R,Rn) mit einer Basis uiC1(R,Rn), ui(0)=ei, i=1,,n.
    Die normierte Fundamentalmatrix ist Y0(t)=(u1,,un).

  • Die Lösungen der inhomogen DGL u(t)=A(t)u(t)+b(t) bilden einen affinen Unterraum u~+VC1(R,Rn) mit einer speziellen Lösung u~. Für die Lösung gilt
    u(t)=Y0(t)u0+0tY0(t)(Y0(s))1b(s)ds (dabei sei t0=0).

  • Ist die DGL autonom, d. h. ist u(t)=Au(t), so gilt Y0(t)=eAt:=n=0tnn!An.

Beispiel: (u1(t)u2(t)) = (u2(t)u1(t)) + (10), d. h. A= (0110). Es gilt A2=E2, d. h.
Y0(t)=n=0tnn!An=k=0t2k+1(2k+1)!A+k=0t2k(2k)!E2=sinh(t)A+cosh(t)E2= (cosh(t)sinh(t)sinh(t)cosh(t)).
Wegen detY0(t)=cosh2(t)sinh2(t)=1 gilt Y01(t)= (cosh(t)sinh(t)sinh(t)cosh(t)) und somit ist die Lösung u(t)=2 (sinh(t)cosh(t)) (01).

Beispiel: Für u(t)=t3u+et, u(0)=1 ist die homogene DGL uh=t3uh, deren Lösung ist uh(t)=et4/4=Y0(t). Die allgemeine Lösungsformel von oben ergibt nun
u(t)=et4/4+et4/40teττ4/4dτ, jedoch kann das Integral analytisch nicht berechnet werden.
Die bestehenden Möglichkeiten sind nun einerseits das Anwenden einer Quadraturformel für das Integral, zum anderen numerische Verfahren für das Ausgangsproblem.

Einzelschrittverfahren

Einzelschrittverfahren:  Angenommen, das Anfangswertproblem besitzt eine eindeutige Lösung uC1(I,Rn). Seien t0:=0 und I:=[0,T] mit T>0.
Ein Schrittweitenvektor ist ein Vektor h:=(h0,,hN1)t[0,T]N mit j=0N1hj=T.
Das Gitter Ih zu h ist Ih:={0=t0,t1,,tN=T} mit tj:=tj1+hj1.
Das Gitter heißt äquidistant, falls h0==hN1. In diesem Fall sei h skalar (h=h0).
Die Gitterweite ist |h|:=maxj=0,,N1hj.
Das Ziel ist die Bestimmung einer Gitterfunktion uh:IhRn. Dabei setzt man uj:=uh(tj) für j=0,,N.

Das Eulersche Polygonzugverfahren

Bemerkung: Zur Vereinfachung setzt man n=1, Ih äquidistant und uh(t0)=u0.
Für die exakte Lösung u des Anfangswertproblems gilt u(t1)=u(t0)+u(t01)(t1t0)
=u0+hf(t01,u(t01)) mit t01[t0,t1] (Taylorformel mit Restglied).
Mittels t01t0 erhält man eine Näherung u1=uh(t1) für u(t1), wobei u1=u0+hf(t0,u0).

explizites Euler-Verfahren: 
Das explizite Euler-Verfahren hat die Iterationsvorschrift uj:=uj1+hf(tj1,uj1).

Beispiel: Für u(t)=t3u+et, u(0)=1 und t[0,1] erhält man schon für geringe N gute Näherungen. Bei u(t)=sin(t)u(t), u(0)=1 (exakte Lösung u(t)=e1cos(t)) und t[0,50] benötigt man schon wesentlich größere Werte für N, um sinnvolle Näherungen zu erzeugen.

Allgemeine Definition, Beispiele

explizites Einschrittverfahren:  Es seien ein Gitter Ih und eine Funktion
ϕC([0,T]2×Rn,Rn) gegeben. Dann heißt das Verfahren uj:=uj1+hj1ϕ(hj1,tj1,uj1), j=1,,N explizites Einschrittverfahren (ESV) und ϕ heißt zugehörige Inkrementfunktion.

Beispiel: Im Euler-Verfahren setzt man u(t01)u(t0)=f(t0,u0).
Man kann dies auch anders approximieren: u(t01)f(t0+h2,u(t0+h2)) mit
u(t0+h2)u(t0)+h2u(t0)=u0+h2f(t0,u0). Daraus ergibt sich die neue Iterationsvorschrift uj:=uj1+hj1f(tj1+hj12,uj1+hj12f(tj1,uj1)), j=1,,N.
Dieses Verfahren nennt sich modifiziertes explizites Euler-Verfahren.

Beispiel: Ein anderes Verfahren ergibt sich wie folgt: u(t1)=u(t0)+(t1t0)u(t01)
=u0+hf(t01,u(t01))=u0+h2(f(t01,u(t01))+f(t01,u(t01)))u0+h2(f(t0,u(t0))+f(t1,u(t1)))
u0+h2(f(t0,u0)+f(t0+h,u0+hf(t0,u0))).
Das sogenannte Verfahren von Heun hat also die Iterationsvorschrift
uj:=uj1+hj12(f(tj1,uj1)+f(tj1+hj1,uj1+hj1f(tj1,uj1))), j=1,,N.

explizites Euler-Verfahren:  Die Inkrementfunktion des expliziten Euler-Verfahrens ist
ϕ(k,t,w):=f(t,w).

modifiziertes explizites Euler-Verfahren:  Die Inkrementfunktion des
modifizierten expliziten Euler-Verfahrens ist ϕ(k,t,w):=f(t+k2,w+k2f(t,w)).

Verfahren von Heun:  Die Inkrementfunktion des Verfahrens von Heun ist
ϕ(k,t,w):=12(f(t,w)+f(t+k,w+kf(t,w))).

Konsistenz, Konvergenz, Stabilität, numerischer Aufwand

globale Fehlerfunktion/globaler Diskretisierungsfehler: 
Die Funktion eh:IhRn mit eh:=u|Ihuh heißt globale Fehlerfunktion.
Der globale Diskretisierungsfehler ist eh¯:=maxj=0,,Neh(tj).

lokale Fehlerfunktion/lokaler Diskretisierungsfehler: 
Die Funktion εh:IhRn mit εh(tj)=1hj(u(tj+1)u(tj)hjϕ(hj,tj,u(tj))) heißt
lokale Fehlerfunktion. Der lokale Diskretisierungsfehler ist εh¯:=maxj=0,,Nεh(tj).

Bemerkung: Der lokale Diskretisierungsfehler gibt den Fehler an, der bei einem Schritt gemacht wird. Er kann als Differenz von der Steigung der exakten Lösung u und der Steigung der Approximation uh interpretiert werden.

Konvergenz:  Das Einzelschrittverfahren heißt konvergent, falls eh¯0 für |h|0.

Konsistenz:  Das Einzelschrittverfahren heißt konsistent zu (AWP), falls εh¯0 für |h|0.

Konsistenzordnung: 
Das Einzelschrittverfahren heißt konsistent zur Ordnung p zu (AWP), falls εh¯=O(|h|p).

numerischer Aufwand:  Der numerische Aufwand ist die Anzahl der Auswertungen von f.

Beispiel: Konsistenz und numerischer Aufwand der bisher betrachteten Verfahren:
explizites Euler-Verfahren: p=1 und 1
modifiziertes Euler-Verfahren: p=2 und 2
Verfahren von Heun: p=2 und 2

Bemerkung: Der Aufwand pro Zeitschritt ist proportional zu p.

Satz (Konsistenz von Einzelschrittverfahren): Seien h[0,T]N ein Schrittweitenvektor, Ih ein Gitter und ϕC([0,T]2×Rn,Rn) die Inkrementfunktion für ein Einzelschrittverfahren (ESV).

  • Das Einzelschrittverfahren ist zu (AWP) konsistent genau dann, wenn
    tIϕ(0,t,u(t))=f(t,u(t)).

  • Seien zusätzlich fCp(I×Rn,Rn) und ϕCp([0,T]2×Rn,Rn).
    Dann ist das Einzelschrittverfahren konsistent mit der Ordnung p zu (AWP) genau dann, wenn tIdidtif(t,u(t))=(i+1)ikiϕ(k,t,u(t))|k=0 für i=0,,p1.

Bemerkung: Was ist der Zusammenhang zwischen dem lokalen Konsistenzfehler εh und dem globalen Fehler eh?

Raum der beschränkten Gitterfunktionen:  Sei Ih ein Gitter zum Schrittweitenvektor h. Die Menge Xh:={vh:Ih{tn=T}Rn|c>0j=0,,N1vh(tj)c} heißt
Raum der beschränkten Gitterfunktionen.
Mit der Norm vh:=maxj=0,,N1vh(tj) ist Xh ein Banachraum isomorph zu RnN.

diskreter Operator:  Seien Ih ein Gitter zum Schrittweitenvektor h und ϕC(I2×Rn,Rn) die Inkrementfunktion für ein Einzelschrittverfahren (ESV). Der Operator Th:XhXh,
(Thvh)(t0):=vh(t0)u0 und (Thvh)(tj):=1hj(vh(tj+1)vh(tj)hjϕ(hj,tj,vh(tj))) für
j=1,,N1 heißt der dem Einzelschrittverfahren (ESV) zugeordnete diskrete Operator.

Bemerkung:
uh ist die Gitterfunktion aus einem Einzelschrittverfahren genau dann, wenn Thuh=0.
Es gilt Th(u|Ih)=O(|h|p), da (Th(u|Ih))(tj)=εh(tj), falls (ESV) kons. mit Ordn. p ist.

Stabilität:  Der Operator Th heißt stabil, falls
c,h¯>0h[0,T]N,|h|<h¯vh(1),vh(2)Xhvh(1)vh(2)cThvh(1)Thvh(2).

Bemerkung: Sei das Einzelschrittverfahren (ESV) stabil. Dann gilt:
Die Lösung uh von Thuh=0 ist eindeutig, denn uhu~hcThuhThu~h=0.
Die Lösung uh von Thuh=0 ist beschränkt, denn
uh=uh0cThuhTh[0]=cc0 für c0:=Thuh.

Satz (Konvergenz von Einzelschrittverfahren I):
Sei ein ESV mit Inkrementfunktion ϕC(I2×Rn,Rn) gegeben. Ist das ESV stabil, so gilt:

  • Ist das ESV konsistent, so ist es auch konvergent.

  • Ist das ESV konsistent zur Ordnung pN, so gilt eh¯=O(|h|p).

Satz (Konvergenz von Einzelschrittverfahren II):
Sei ein ESV mit Inkrementfunktion ϕC(I2×Rn,Rn) gegeben. Außerdem existiere für h¯ fest eine Konstante M>0 mit k(0,h¯)tIw,w~Rnϕ(k,t,w)ϕ(k,t,w~)Mww~
(globale Lipschitz-Bedingung an ϕ im dritten Argument).

  • Das ESV ist stabil.

  • Ist das ESV konsistent, so ist es auch konvergent.

  • Ist das ESV konsistent zur Ordnung p, so existiert eine Konstante c>0, sodass für alle Gitter Ih mit |h|<h¯ die Abschätzung eh¯ccs|h|p gilt, wobei cs:=eMT(T+1) die Stabilitätskonstante und I=[0,T] ist.

Bemerkung: Die Abschätzung von (iii) ist bzgl. der Stabilitätskonstanten cs bestmöglich, d. h. auf I=[0,) ist nicht mit gleichmäßiger Konvergenz zu rechnen.

Beispiel: Als Beispiel betrachtet man das AWP u(t)=au(t) mit u(0)=1 und a>0. Für die Lösung u(t)=eat ergibt sich bei Anwendung des expliziten Euler-Verfahrens mit äquidistantem Gitter uj=(1+ah)j1=(1+ah)tj/h1 (mit tj=jh), also eh(tj)=eatj(1+ah)tj/h1.

Satz (Konvergenz von Einzelschrittverfahren III):
Sei ein ESV mit Inkrementfunktion ϕC(I2×Rn,Rn) gegeben. Außerdem existiere für h¯ und ε>0 fest eine Konstante M>0 mit
k(0,h¯)tIw,w~{vRn|tIvu(t)ε}ϕ(k,t,w)ϕ(k,t,w~)Mww~
(lokale Lipschitz-Bedingung an ϕ im dritten Argument).
Dann gelten (i), (ii) und (iii) aus obigem Satz:

  • Das ESV ist stabil.

  • Ist das ESV konsistent, so ist es auch konvergent.

  • Ist das ESV konsistent zur Ordnung p, so ist es auch konvergent zur Ordnung p.

Explizite Runge-Kutta-Verfahren

Bemerkung: Seien pN0 und ein Anfangswertproblem (AWP) mit einer Lösung
uCp+1(I,Rn) vorgegeben (dies ist z. B. der Fall für fCp(I,Rn)).
Kann man nun systematisch ein Einzelschrittverfahren mit Konsistenzordnung p konstruieren?

Beispiel: Das Heun-Verfahren uj+1=uj+hj2(f(tj,uj)+f(tj+hj,uj+hjf(tj,uj))) mit Inkrementfunktion ϕ(k,t,w)=12(f(t,w)+f(t+k,w+kf(t,w))) erreicht durch iterative Auswertung von f eine höhere Konsistenzordnung (nämlich 2). Es gehört zu den einfachsten expliziten Runge-Kutta-Verfahren.

explizites Runge-Kutta-Verfahren:  Seien rN, α2,,αrR, γ1,,γrR und βij für i=2,,r, j=1,,r1 und i>j gegeben.
Das Einzelschrittverfahren (ESV) mit ϕ(k,t,w):=i=1rγiKi(k,t,w) und
K1(k,t,w):=f(t,w),
K2(k,t,w):=f(t+α2k,w+kβ21K1(k,t,w)),
…,
Kr(k,t,w):=f(t+αrk,w+ks=1r1βrsKs(k,t,w))
heißt allgemeines explizites Runge-Kutta-Verfahren der Stufe r.

Butcher-Tableau:  Die Koeffizienten eines allgemeinen Runge-Kutta-Verfahrens können in der Form einer Tabelle (Butcher-Tableau) zusammengefasst werden:

α2β21αrβr1βr,r1γ1γr1γr

Beispiel: Das explizite Euler-Verfahren ϕ(k,t,w)=f(t,w) (Stufe 1), das modifizierte Euler-Verfahren ϕ(k,t,w)=f(t+k2,w+k2f(t,w)) (Stufe 2) und das Verfahren von Heun
ϕ(k,t,w)=12(f(t,w)+f(t+k,w+kf(t,w))) (Stufe 2) besitzen folgende Butcher-Tableaus:

1121201111212

Bemerkung: Setzt man α1:=0, so kann man die Koeffizientenfunktionen Ki iterativ bestimmen durch die Formel Ki(k,t,w)=f(t+αik,w+kj=1i1βijKj(k,t,w)) für i=1,,r.

Beispiel: Für K1=f(tj,uj), K2=f(tj+hj2,uj+12hjK1), K3=f(tj+hj2,uj+12hjK2), K4=f(tj+hj,uj+hjK3) ergibt sich ein Runge-Kutta-Verfahren mit der Inkrementfunktion uj+1=uj+hj6(K1+2K2+2K3+K4). Es heißt klassisches Runge-Kutta-Verfahren und besitzt folgendes Butcher-Tableau:

121212012100116131316

Bemerkung: Man kann sich fragen, wieviele Runge-Kutta-Verfahren r-ter Ordnung auch eine Konsistenzordnung von r besitzen.
Für den Fall r=2 ergibt obiger Satz (Konsistenz von Einzelschrittverfahren) die Bedingungen f(t,w)=ϕ(0,t,w) und ddtf(t,u(t))=2kϕ(k,t,u(t))|k=0.
Es gilt ϕ(0,t,w)=γ1f(t,w)+γ2f(t,w) und ddtf(t,u(t))=ft+fududt=ft+fuf(t,u(t)).
Für die Ableitung von ϕ gilt ϕ(k,t,w)=γ1f(t,w)+γ2f(t+α2k,w+kβ21f(t,w)), also
2ϕk|k=0=2γ2(α2ft+β21f(t,u(t))fu).
Aus Koeffizientenvergleich ergibt sich das nicht-lineare Gleichungssystem 1=γ1+γ2, 2γ2α2=1, 2γ2β21=1. Für γ20 kann man die drei Gleichungen mit vier Unbekannten mit γ2 als Parameter auflösen und erhält γ1=1γ2, α2=12γ2 und β21=12γ2. Das Butcher-Tableau lautet

12γ212γ21γ2γ2. Für γ2=12 erhält man das Heun-Verfahren und für γ2=1 das modifizierte Euler-Verfahren.

Bemerkung: Allgemein muss man ein nicht-lineares Gleichungssystem lösen. Die Konsistenzordnung eines r-stufigen Runge-Kutta-Verfahrens ist nach oben durch r beschränkt.
Leider gilt i. A. nicht, dass r die maximal erreichbare Konsistenzordnung ist:

r123456789pmax(r)12344566r2

Bemerkung: Ein Runge-Kutta-Verfahren der Stufe r ist konsistent genau dann, wenn
i=1rγi=1 gilt, denn aufgrund Ki(0,t,w)=f(t,w) für i=1,,r gilt
ϕ(0,t,w)=i=1rγiKi(0,t,w)=f(t,w)i=1rγi=!f(t,w).

Implizite Runge-Kutta-Verfahren

Bemerkung: Man spricht von einem impliziten Einzelschrittverfahren, falls die Inkrementfunktion ϕ auch von ui+1=uh(ti+1) abhängt, d. h. ui+1=ui+hiϕ(hiti,ui,ui+1).
Die Vorteile sind die verbesserte Stabilität und eine höhere mögliche Konsistenzordnung von bis zu 2r. Der Nachteil ist natürlich der höhere numerische Aufwand, da man pro Zeitschritt ein in der Regel nicht-lineares Gleichungssystem lösen muss.

Beispiel: Das implizite Euler-Verfahren ist gegeben durch ui+1=ui+hif(ti+1,ui+1).

implizites Runge-Kutta-Verfahren: 
Seien rN, α1,,αrR, γ1,,γrR und bij für i,j=1,,r gegeben.
Das Einzelschrittverfahren (ESV) mit ϕ(k,t,w):=i=1rγiKi(k,t,w) heißt allgemeines
implizites Runge-Kutta-Verfahren
der Stufe r, falls das nicht-lineare Gleichungssystem
K1(k,t,w):=f(t+α1k,w+ks=1rb1sKs(k,t,w)),
…,
Kr(k,t,w):=f(t+αrk,w+ks=1rbrsKs(k,t,w)),
erfüllt ist. Die Koeffizienten können analog zum expliziten Fall in einem Butcher-Tableau zusammengefasst werden:

α1b11b1rαrbr1brrγ1γr

Beispiel: Beim impliziten Euler-Verfahren uj+1=uj+hjf(tj+1,uj+1) ist K1=f(tj+1,uj+1), d. h. K1=f(tj+1,uj+hjK1). In der Regel ist pro Zeitschritt eine nicht-lineare Gleichung zu lösen. Man kann z. B. eine einfache Iteration K1(+1)=f(tj+1,uj+hhK1()) bzw.
uj+1(+1)=uj+hjf(tj+1,uj+1()) lösen oder das Newton-Verfahren anwenden.
Das Butcher-Tableau ist folgendes:

111

Beispiel: Allgemeine Runge-Kutta-Verfahren der Stufe r=1 besitzen im skalaren Fall n=1 die nicht-lineare Gleichung K1=f(t+α1k,w+kb11K1). Man nimmt an, dass die Gleichung eindeutig lösbar ist mit K1C1(I2×R,R).
Für die Konsistenz muss ϕ(0,t,w)=γ1K1(0,t,w)=f(t,w) gelten, das stimmt für γ1=1 (für α1=b11=1 erhält man das implizite Euler-Verfahren). Differentiation von obiger Gleichung in k=0 ergibt kϕ(0,t,w)=ft(t,w)α1+fw(t,w)b11K1(0,t,w) und
ddtf(t,u(t))=ft(t,u(t))+fw(t,u(t))f(t,u(t)). Nach dem Konsistenzsatz muss für p=2 gelten, dass 2kϕ(0,t,u(t))=ddtf(t,u(t)), also α1=b11=12.
Konkret erhält man also uj+1=uj+hjK1=uj+hjf(tj+12hj,uj+12hjK1)
=uj+hjf(12(tj+tj+1,uj+12(uj+1uj))=uj+hjf(12(tj+tj+1),12(uj+uj+1)),
da K1=1hj(uj+1uj).

Beispiel: implizites Runge-Kutta-Verfahren der Stufe r=2 und Ordnung p=4:

(33)/61/41/43/6(3+3)/61/4+3/61/41/21/2

Bemerkung: Jedes implizites ESV lässt sich für |h| hinreichend klein als explizites Verfahren darstellen.

Beispiel: Wendet man das implizite Euler-Verfahren auf u=au an, so erhält man
ui+1=ui+ahui+1, also ui+1=11ahui=ui+(11ah1)ui=ui+ha1ahui.

Beispiel: Für die halbimpliziten Runge-Kutta-Verfahren gilt bis=0 für i<s, also
K1(k,t,w):=f(t+α1k,w+kb11K1(k,t,w)),
…,
Kr(k,t,w):=f(t+αrk,w+ks=1rbrsKs(k,t,w)),
d. h. die einzelnen Gleichungen sind nacheinander lösbar. Das Butcher-Tableau hat dann folgende Form:

α1b110αrbr1brrγ1γr

Zusammenhang zwischen Runge-Kutta-Verfahren und Quadraturformeln

Bemerkung: Der Zusammenhang zwischen Runge-Kutta-Verfahren und Quadraturformeln gibt einen weiteren Weg zur systematischen Konstruktion von Runge-Kutta-Verfahren zu einer vorgegebenen Konsistenzordnung p.

Bemerkung: Gegeben sei das allgemeine Runge-Kutta-Verfahren uj+1=uj+hji=1rγiKi,
Ki=f(tj+αihj,uj+hjs=1rbisKs).
Im Folgenden wird versucht, eine notwendige Bedingung für die Konsistenzordnung p herzuleiten. Betrachtet man das Anfangswertproblem u(t)=g(t), tI, u(0)=u0Rn mit g:IRn (Lösung u(t)=u0+t0tg(τ)dτ), so ergibt das Runge-Kutta-Verfahren
1hj(uj+1uj)=i=1rγig(tj+αihj). Der Konsistenzfehler ist laut Definition
εh(tj)=1hj(u(tj+1)u(tj)hji=1rγig(tj+αihj))=1hj(tjtj+1g(t)dthji=1rγig(tj+αihj)), d. h. um die Konsistenzordnung p zu erreichen, muss
|tjtj+1g(t)dthji=1rγig(tj+αihj)|=O(hjp+1) gelten.
Dies ist ein Quadraturproblem (Gewichte γi, Stützstellen tj+αihj). Damit können die Koeffizienten α1,,αr und γ1,,γr bestimmt werden.

Bemerkung: Das Einsetzen der exakten Lösung in das Runge-Kutta-Verfahren ergibt
u(tj+1)u(tj)hji=1rγiKi(hj,tj,u(tj)). Daraus folgt tjtj+1u(t)dthji=1rγiKi(hj,tj,u(tj)).
Wegen der Gleichung von eben sollte für schon gegebene αi und γi gelten, dass Kiu(tj+αihj) für i=1,,r. Aus der Definition der Ki folgt damit Ki=f(tj+αihj,uj+hjs=1rbisKs)u(tj+αihj)=f(tj+αihj,u(tj+αihj)). Daraus folgt u(tj+αihj)u(tj)+h(tj)s=1rbisKs, d. h. tjtj+αihju(t)dts=1rbisu(tj+αshj) für i=1,,r. Somit erhält man ein Quadraturproblem, mit dem sich die bis bestimmen lassen (i,s=1,,r).
Dies motiviert den folgenden Satz.

Satz (Butcher, Kunzmann, 1969): Es sei ein Runge-Kutta-Verfahren mit αi,γiR für i=1,,r und bijR für i,j=1,,r gegeben. Für p,qN seien die Koeffizienten so gewählt, dass für alle g1Cp+1(I,Rn) und g2Cq+1(I,Rn) gilt

  • |1hjtjtj+1g1(t)dts=1rγsg1(tj+αshj)|=O(hjp) für j=0,,N1 und

  • |1hjtjtj+αihjg2(t)dts=1rβisg2(tj+αshj)|=O(hjq) für j=0,,N1 und i=1,,r.

Dann ist das Runge-Kutta-Verfahren konsistent mit der Ordnung min{p,q+1}.

Exaktheit einer Quadraturformel:  Es seien gC([0,1],R) und τ(0,1] gegeben.
Sei Q(g):=i=1rγig(αi) eine Quadraturformel für das Integral 0τg(t)dt, wobei αi[0,1] und γiR für i=1,,r. Q heißt vom Grad exakt, falls Q(p)0τp(t)dt=0 für alle pP
(P Menge der Polynome vom Grad ).

Satz (Fehler einer Quadraturformel mit Peano-Kern): Seien N, gC+1([0,1],R),
τ(0,1] und Q eine Quadraturformel, die vom Grad exakt ist.
Dann gilt 0τg(t)dt=Q(g)+01π+1(t)g(+1)(t)dt, wobei π+1 der Peano-Kern
π+1(t):=1(+1)!(((τt)+)+1(+1)i=1rγi((αit)+)) ist mit
t[0,1] und α+(t)=max{α(t),0}.

Legendre-Polynom:  Für mN0 ist das Legendre-Polynom pm vom Grad m gegeben durch pm(t):=m!(2m)!dmdtm(t21)m für tR.

Beispiel: Es gilt p0(t)=1, p1(t)=t, p2(t)=t213 usw.

Lemma (Nullstellen und Orthogonalität der Legendre-Polynome):

  • Das Legendre-Polynom pm besitzt paarweise verschiedene Nullstellen
    ϱ1,,ϱm mit 1<ϱ1<<ϱm<1.

  • Für m,nN mit mn gilt 11pm(t)pn(t)dt=0.

Satz (Gauß-Quadratur): Seien gC([1,1],R) und Q(g):=i=1mωig(ϱi) die
Gauß-Quadraturformel mit den Stützstellen ϱi (Nullstellen des Legendre-Polynoms pm) und den Gewichten ωi:=11(j=1,jimtϱjϱiϱj)dt (Integrale für Lagrange-Polynome) für mN.
Dann gilt Q(p)=11p(t)dt für alle pP2m1,
d. h. die Gauß-Quadratur ist exakt vom Grad 2m1.

Bemerkung: Nun kann man analysieren, wie gut das Runge-Kutta-Verfahren ist, das durch die Gauß-Quadratur bestimmt wird. Dazu wendet man den Satz von Butcher und Kunzmann an.

  • Mit t=tj+hjτ, τ[0,1] gilt 1hjtjtj+1g1(t)dt=01g1(tj+hjτ)dτ
    =i=1rω~ig1(tj+hjϱ~i)+01π2r(τ)g1(2r)(tj+hjτ)dτ, da die Gauß-Quadratur exakt vom Grad 2r1 ist. Daraus folgt |1hjtjtj+1g1(t)dti=1rω~ig1(tj+hjϱ~i)|
    maxτ[0,1]|π2r(τ)|01|g1(2r)(tj+hjτ)|dτc|h|2r aufgrund der Beschränktheit von π (bei jeder Ableitung von g1 kommt ein Faktor hj hinzu). Dabei ist γi:=ω~i=ωi2 und αi:=ϱ~i=ϱi+12 für i=1,,r, weil 01g(τ)dτ=1211g(z+12)dz12i=1rωig(ϱi+12).

  • Analog wie eben ist 1hjtjtj+αihjg2(t)dt=0αig2(tj+hjτ)dτ
    =s=1rω^sg2(tj+hjϱ^s)+01π2r(τ)g2(2r)(tj+hjϱ^s)dτ. Daraus folgt wieder
    |1hjtjtj+αihjg2(t)dts=1rω^sg2(tj+hjϱ^s)|c|h|2r mit
    βis:=ω^s=αiωs2 und ϱ~s=αiϱs+12 wegen 0αg(τ)dτ=α211g(αz+12)dz.

Somit ergibt sich eine Konsistenzordnung von p=min{2r,2r+1}=2r.

Mehrschrittverfahren

Definitionen und Beispiele

Bemerkung: Um die Genauigkeit von Einzelschrittverfahren zu erhöhen, verwendet man nicht nur die letzte, sondern die letzten k Approximationen.

Mehrschrittverfahren:  Seien ψC(Ik+2×Rn(k+1),Rn) und kN.
Weiter seien a0,,akR und u0=u(t0),u1,,uk1Rn gegeben. Das Verfahren
1h(a0uj+a1uj+1++akuj+k)=ψ(h,tj,,tj+k,uj,,uj+k) mit j=0,,Nk heißt
k-Mehrschrittverfahren (k-MSV) mit Verfahrensfunktion ψ. (Das Gitter Ih ist also äquidistant.)

Bemerkung: Falls ψ nicht von uj+k abhängt und ak0 gilt, so heißt das k-MSV explizit.
Ein explizites 1-MSV ist ein explizites Einzelschrittverfahren.

Bemerkung: Um die Verfahrensgleichung lösen zu können, müssen zunächst die Startwerte u0,,uk1 bekannt sein. Diese sollte mit einem ESV derselben Konsistenzordnung berechnet werden.

lineares MSV:  Falls die Verfahrensfunktion ψ von der Form ψ(h,tj,,tj+k,uj,,uj+k)=i=0kbif(tj+i,uj+i) mit b0,,bkR ist, so heißt das zugehörige k-MSV linear.
Lineare MSV haben also die Form 1hi=0kaiuj+i=i=0kbif(tj+i,uj+i)

Bemerkung: Auch MSV lassen sich durch Quadraturformeln herleiten. Die Integralgleichung u(t)=u(s)+stf(r,u(r))dr mit t>s und t,sI ist äquivalent zu u(t)=f(t,u(t)), speziell gilt u(tj+k)=u(tj+k1)+tj+k1tj+kf(r,u(r))dr=u(tj+k1)+tj+k1tj+ku(r)dr.

Beispiel: Verwendet man die Trapezregel
u(tj+1)=u(tj)+tjtj+1f(r,u(r))dru(tj)+12h(f(tj,u(tj))+f(tj+1,u(tj+1))) (also k=1), so ergibt sich das Trapezverfahren uj+1:=uj+12h(f(tj,uj)+f(tj+1,uj+1)).

Bemerkung: Eine Idee für weitere Verfahren ist eine bessere Approximation der Integralgleichung durch Ersetzung des Integranden u(r) durch ein Interpolationspolynom.
Die Interpolationspolynome pjPk1, j=0,,Nk sind eindeutig bestimmt durch die k Bedingungen pj(tj+i):=u(tj+i)=f(tj+i,u(tj+i)) für i=0,,k1. Man erhält die veränderte Integralgleichung u(tj+k)u(tj+k1)+tj+k1tj+kpj(r)dr.
Da man dafür allerdings die exakte Lösung u benötigt, kann man auch p~jPk1 verwenden, die analog definiert sind durch p~j(tj+i):=f(tj+i,uj+i) für i=0,,k1. Verwendet man p~j statt pj in der Integralgleichung, so erhält man ein explizites k-MSV
1h(uj+kuj+k1)=ψ(h,tj,,tj+k,uj,,uj+k1):=1htj+k1tj+kp~j(r)dr.

Beispiel: Ein Beispiel für ein so erhaltenes lineares 4-Mehrschrittverfahren mit p=4 ist
1h(uj+4uj+3)=124(55f(tj,uj)59f(tj+1,uj+1)+37f(tj+2,uj+2)9f(tj+3,uj+3)).
Es heißt Adams-Bashforth-Verfahren der Stufe k=4.

Bemerkung: Bei jedem Zeitschritt ist nur eine neue Auswertung von f notwendig
(und zwar in (tj+k1,uj+k1)).

Beispiel: Ein implizites Verfahren lässt sich analog konstruieren, nur bezieht man dabei
tj+k,uj+k als Stützpunkte für die Interpolation ein.
Diese Verfahren heißen Adams-Moulton-Verfahren.
Ein Beispiel für k=4 und p=5 ist 1h(uj+4uj+3)
=1720(251f(tj+4,uj+4)+646f(tj+3,uj+3)269f(tj+2,uj+2)+106f(tj+1,uj+1)19f(tj,uj).

Bemerkung: Bei den sog. Prädiktor-Korrektor-Verfahren kombiniert man implizite und explizite Verfahren. Seien also 1hi=0kaiuj+i=ψ1(h,tj,,tj+k,uj,,uj+k1) ein explizites und 1hi=0kαiuj+i=ψ2(h,tj,,tj+k,uj,,uj+k) ein implizites k-MSV.
Man berechnen nun zuerst den Prädiktor uj+k(p) mit dem expliziten MSV, d. h.
1hi=0k1aiuj+i+1hakuj+k(p)=ψ1(h,tj,,tj+k,uj,,uj+k1). Anschließend berechnet man uj+k mit dem impliziten Verfahren, also 1hi=0kαiuj+i=ψ2(h,tj,,tj+k,uj,,uj+k1,uj+k(p)).
Man muss also keine nicht-linearen Gleichungen lösen, sondern man verwendet den Prädiktor als Schätzwert für den wahren Wert uj+k.
Alternativ lässt sich der Prädiktor auch als Startwert für eine Fixpunktiteration verwenden, d. h. uj+k(0):=uj+k(p) und 1hi=0k1αiuj+i+1hαkuj+k(m+1):=ψ2(h,tj,,tj+k,uj,,uj+k1,uj+k(m)).

Konsistenz und Konvergenz von Mehrschrittverfahren

Fehler von linearen Mehrschrittverfahren: 
Es sei uh:IhRn durch ein lineares k-MSV gegeben.
eh:=u|Ihuh ist die globale Fehlerfunktion.
eh¯:=maxj=0,,Neh(tj) ist der globale Diskretisierungsfehler.
εh(tj+k):=1hi=0kaiu(tj+i)i=0kbif(tj+i,u(tj+i)), j=0,,Nk, ist die lokale Fehlerfunkt.
εh¯:=maxj=0,,Nkεh(tj+k) ist der lokale Diskretisierungsfehler.

Bemerkung: Die Koeffizienten a0,,ak,b0,,bk sollten so bestimmt werden, dass
εh¯=O(hp). Dafür betrachtet man εh(tj+k)=1hi=0kaiu(tj+ih)i=0kbiu(tj+ih) und setzt für p(i):=1!i die Taylor-Entwicklungen u(tj+ih)==0php(i)u()(tj)+O(hp+1) bzw.
u(tj+ih)==0php(i)u(+1)(tj)+O(hp+1)==1ph1p1(i)u()(tj)+O(hp)
==1ph1p(i)u()(tj)+O(hp) ein. Daraus folgt dann
εh(tj+k)=1hi=0kai(=0php(i)u()(tj))i=0kbi(=1ph1p(i)u()(tj))+O(hp)
==0ph1u()(tj)(i=0kaip(i)i=0kbip(i))+O(hp).
Verschwindet der Ausdruck in Klammern, so hat das Verfahren die Konsistenzordnung p. Das beweist folgenden Satz.

Satz (Konsistenz von MSV): Falls die Koeffizienten eines linearen k-MSV
a0,,ak,b0,,bkR die Bedingungen i=0kaip(i)=i=0kbip(i) für =0,,p erfüllen, so besitzt das MSV die Konsistenzordnung p.
Dabei ist p(i):=1!i und p(i):=p1(i)=1(1)!i1 für 1 bzw. p0(i):=0.

Bemerkung: Diese Bedingungen entsprechen einem LGS mit p+1 Gleichungen und 2(k+1) Unbekannten. Da die Lösung a0==ak=b0==bk=0 keinen Sinn ergibt, ergänzt man manchmal die Normierungsbedingung i=0kbi=1.
Damit das Gleichungssystem nicht überbestimmt ist, soll es höchstens so viele Gleichungen wie Variablen geben. Mit der Normierungsbedingung ist dann p+22(k+1), d. h. die Konsistenzordnung p ist durch 2k nach oben beschränkt.
Bei expliziten Verfahren ist bk=0, d. h. es gibt eine Variable weniger. Hier ist p+22k+1, also ist die Konsistenzordnung p durch 2k1 nach oben beschränkt.

Beispiel: Für k=1 soll p=2 erreicht werden, d. h. die Gleichungen a0+a1=0, a1=b0+b1, 12a1=b1 und b0+b1=1 sollen erfüllt werden. Daraus folgt a0=1, a1=1, b0=12 und b1=12. Man erhält also die Trapezregel 1h(uj+uj+1)=12f(tj,uj)+12f(tj+1,uj+1).

Stabilität von Mehrschrittverfahren

erzeugende Polynome: 
Sei ein lineares k-Mehrschrittverfahren 1hi=0kaiuj+i=i=0kbif(tj+i,uj+i) gegeben.
Dann heißen die Polynome ϱ(z):=i=0kaizi und σ(z):=i=0kbizi erzeugende Polynome des MSV (zC).

alternative Schreibweise von linearen MSV:  Sei E der Vorwärts-Shift-Operator, d. h. Eyj:=yj+1. Dann lässt sich das lineare k-Mehrschrittverfahren
1hi=0kaiuj+i=i=0kbif(tj+i,uj+i) auch durch die erzeugenden Polynome in der Form
1hϱ(E)uj=σ(E)fj mit fj=f(tj,uj) schreiben, wobei p(E)yj:=i=0kpiyj+i mit einem Polynom p(z)=i=0kpizi.

Bemerkung: Nicht jedes konsistente lineare MSV ist konvergent. Es wird eine zusätzliche Stabilitätsbedingung benötigt.

Beispiel: Ein Beispiel für ein instabiles lineares 2-MSV mit Konsistenzordnung p=3 ist
1h(ui+2+4ui+15ui)=4fi+1+2fi. Die erzeugenden Polynome sind dabei ϱ(z)=z2+4z5 und σ(z)=4z2. Man wendet das MSV auf das triviale Anfangswertproblem u=0, u(0)=1 (d. h. die Lösung ist u(t)1) an.
Sei u1=1+εh leicht gestört. Daraus ergibt sich die Drei-Term-Rekursion ui+2+4ui+15ui=0 (rechte Seite verschwindet wegen f0) mit den Startwerten u0=1 und u1=1+εh.
Für spezielle Lösungen betrachtet man die Nullstellen z1=1 und z2=5 des erzeugenden Polynoms ϱ(z). Setzt man ui=z1i an, so ist z1i+2+4z1i+15z1i=0 genau dann, wenn z1iϱ(z1)=0. Wegen ϱ(z1)=0 ist ui=z1i eine spezielle Lösung der Rekursion, analog ui=z2i.
Für die allgemeine Lösung setzt man ui=Az1i+Bz2i an, also ui+2+4ui+15ui=0 genau dann, wenn Az1iϱ(z1)+Bz2iϱ(z2)=0. Die Parameter A und B ergeben sich aus den Startbedingungen 1=u0=Az10+Bz20=A+B und 1+εh=u1=Az11+Bz21=A5B.
Daraus ergibt sich A=1+εh6 und B=εh6. Somit ist die allgemeine Lösung der Rekursion
ui=Az1i+Bz2i=1+εh6εh6(5)i. Für den Fall ε=0 kommt die exakte Lösung heraus. Ist allerdings u1 leicht gestört (ε>0), so wird der Fehler durch den Faktor (5)i verstärkt, d. h. das MSV ist instabil.

Bemerkung: Die Vorgehensweise lässt sich auf allgemeine k-MSV verallgemeinern. Durch Anwendung der Testgleichung u=0, u(0)=1 erhält man die homogene Rekursion bzw. Differenzengleichung a0uj++akuj+k=0 für j=0,,Nk mit Startwerten u0,,uk1.

Satz (Lösungen der homogenen Rekursion): Sei λC eine m-fache Nullstelle des erzeugenden Polynoms ϱ(z), d. h. ϱ(λ)=ϱ(λ)==ϱ(m1)(λ)=0. Dann gilt:

  • ui(1):=λi,  ui(2):=iλi1,  …,  ui(m):=Dm1λi=i(i1)(im+2)λim+1
    sind spezielle Lösungen der homogenen Rekursion.

  • Die allgemeine Lösung der homogenen Rekursion ist eine Linearkombination der insgesamt k speziellen Lösungen aus (i).
    (Für jede Nullstelle λ von ϱ(z) erhält man entsprechend der Vielfachheit viele spezielle Lösungen, d. h. insgesamt Grad(ϱ)=k viele Lösungen.)

Bemerkung: Sei λ eine Nullstelle von ϱ(z). Dann gilt für |λ|>1, dass {ui}={λi} exponentiell wächst, und für |λ|<1, dass {ui}={λi} exponentiell fällt.
Für |λ|=1 und Vielfachheit von λ ist |ui(1)|=|λ|i=1 und ui()=i(i1)(i+2)λi+1 wächst polynomial für 2.

stabil:  Ein k-Mehrschrittverfahren heißt stabil, falls alle Nullstellen des Polynoms ϱ(z) im abgeschlossenen Einheitskreis liegen und diejenigen auf dem Rand nur einfach sind, d. h.
ϱ(λ)=0|λ|1 und (ϱ(λ)=0|λ|=1)ϱ(λ)0.
Wegen der Testgleichung u=0 spricht man auch von Nullstabilität oder D-Stabilität (nach Dahlquist).

stark/schwach stabil:  Das k-MSV heißt stark stabil, falls für alle Nullstellen außer λ=1 gilt, dass |λ|<1. Ansonsten heißt das k-MSV schwach stabil.

Bemerkung:
Bei konsistenten k-MSV ist λ=1 immer eine Nullstelle von ϱ(z), denn ϱ(1)=i=0kai=0.
Die Adams-Verfahren (Adams-Bashforth und Adams-Moulton) sind stark stabil, denn hier ist ak=1, ak1=1 und ak2==a0=0, d. h. ϱ(z)=zkzk1=zk1(z1).
λ=1 ist einfache Nullstelle, während λ=0 eine (k1)-fache Nullstelle ist.

Satz (Dahlquist-Barriere – maximale Konvergenzordnung stabiler linearer MSV):
Ein lineares k-Mehrschrittverfahren 1hi=0kaiuj+i=i=0kbif(tj+i,uj+i), das obige Stabilitätsbedingung erfüllt, hat maximal die Konvergenzordnung

  • k+2 für k gerade,

  • k+1 für k ungerade und

  • k für bkak0 (insbesondere für explizite Verfahren).

Die Ordnung k+2 kann nur erzielt werden, wenn alle Nullstellen von ϱ(z) auf dem Rand des Einheitskreises liegen.

Beispiel: Für k=1 wird die maximale Konvergenzordnung p=2 von der Trapezformel erreicht. Für k=2 wird die maximale Konvergenzordnung p=4 vom Milne-Simpson-Verfahren ui+1=ui1+h3(fi1+4fi+fi+1) erreicht (schwach stabil).
Das Adams-Bashforth-Verfahren (k=4) ist explizit und erreicht daher nur die Konvergenzordnung p=4. Das Adams-Moulton-Verfahren ist stark stabil (kann nicht Ordnung k+2 erreichen) und erreicht die Konvergenzordnung 5=k+1.

Satz (Konvergenz von MSV): Falls ein lineares MSV die Konsistenzordnung p hat und obige Stabilitätsbedingung erfüllt, so ist es auch konvergent mit Ordnung p.

Adaptive Schrittweitensteuerung

Bemerkung: Sei ein Einzelschrittverfahren zur Lösung des Anfangswertproblems (AWP) gegeben. Für ein gegebenes Gitter Ih sei T(Ih) der numerische Aufwand zur Lösung des ESV auf Ih (Rechenzeit). Außerdem sei TOL eine gegebene Fehlertoleranz.
Die Aufgabe ist nun, ein Gitter Ihopt zu finden mit eh¯TOL und T(Ihopt)T(Ih) für alle Gitter Ih mit eh¯TOL.
Man weiß nicht, ob Ihopt überhaupt existiert oder ob es eindeutig ist. Durch die sog. adaptive Schrittweitensteuerung versucht man, eine möglichst gute Approximation von Ihopt zu finden.

Satz (Fehlerentwicklung): Seien uCp+2(I,Rn) eine Lösung von (AWP) und Ih ein Gitter. Außerdem sei ein stabiles ESV mit Inkrementfunktion ϕCp+1(I2×Rn,Rn), uh(0)=u0 und Konsistenzordnung p gegeben.
Dann existiert eine Funktion e0C2(I,Rn) mit e0(0)=0 und
uh(uhpe0)|Ih=O(hp+1).
Es gibt zusätzlich eine Funktion e1C3(I,Rn) mit e1(0)=0 und
uh(uhpe0hp+1e1)|Ih=O(hp+2).

Bemerkung: Führt man für ein festes tIh das Hilfsproblem v(s)=f(s,v(s)) für s>t und v(t)=uh(t) ein, so gilt uh(t+h)v(t+h)=hpe0(t+h)+hp+1e1(t+h)+O(hp+2)=hp(e0(t)+he0(t))+hp+1(e1(t)+he1(t))+O(hp+2)=hp+1e0(t)+O(hp+2) wegen e0(t)=e1(t)=0. Analog gilt uh/2(t+h)v(t+h)=(h2)phe0(t)+O(hp+2).
Somit ist uh(t+h)uh/2(t+h)=(hp+1(h2)ph)e0(t)+O(hp+2),
d. h. he0(t)=12p1(h2)p(uh(t+h)uh/2(t+h))+O(h2).
Man erhält also für den Fehler der halben Gitterweite die Formel
uh/2(t+h)v(t+h)=Δh+O(hp+2) mit Δh:=12p1(uh(t+h)uh/2(t+h)) einem Fehlerschätzer, der nur aus berechenbaren Größen besteht. Man definiert nun den relativen Fehlerschätzer Δ~h:=Δhmax{1,uh} und kann daraus einen selbstadaptiven Algorithmus erstellen.

selbstadaptiver Algorithmus mit (h,h/2)-Gittersteuerung: 
Startschrittweite h0[0,T],  minimale und maximale Schrittweite hmin<h0<hmax,
Fehlertoleranz TOL>0,  Verkleinerungs-/Vergrößerungsfaktoren kmin<1 und kmax>1,
Verfahrensordnung pN

t:=0;u:=u(0);h:=h0;while(t<T){|Δ~h|:=TOL+1;while(|Δ~h|>TOL){v:=u+hϕ(h,t,u);z:=u+h/2ϕ(h/2,t,u);w:=z+h/2ϕ(h/2,t+h/2,z);|Δ~h|:=12p1|vw|max{1,u};h:=max{hmin,kminh};ifh=hmin{return;}}u:=w;h:=min{hmax,kmaxh};t:=t+h;}