Langzeitverhalten von ODE – L¨osungen

(1)

Mathematik f¨ur Informatiker III Euler Verfahren f¨ur Systeme von ODEs

Bemerkung zum Langzeitverhalten

Häufig ist von Interesse (z.B. in der Klimavorhersage), wie sich Lösungen y(t) der ODE ˙y=F(y) fürsehr grosse tqualitativverhalten, und zwar unabhängig vom Anfangswerty(t0) =y0.

D.h. man will wissen, ob das dynamische System sich einschwingt, einen Gleichgewichtszutand erreicht, zuf¨alliges (d.h. chaotisches) Verhalten o.¨a.

zeigt.

Im folgenden machen wir Aussagen f¨ur autonome Systeme der Zustandsraumdimensionn, die entspechend auch f¨ur nichtautonome Systeme der Dimensionn−1 gelten.

– 79–

(I) Falls

n

= 1 muss und sonst (

n>

1) kann einer der beiden folgenden F¨alle eintreten:

(a)y(t) strebt einem station¨aren Grenzwerty_∞= lim

t→∞y(t) zu Beispiel: ˙y=λ(y−a), a∈R, λ <0,y0beliebig

t y

y_∞

y(t) =c e^λt+a,c<0 y(t) =c e^λt+a,c>0

– 80–

(b)y(t) explodiert (blow up)

tlim→t∗ky(t)k=∞ f¨ur endliche Zeitt^∗(kritische Zeit) Beispiel: y˙=y² mit y(0) =y0>0

=⇒ dy y²=dt =⇒

Z 1 y²dy=

Z

dt =⇒ −1

y=t+c =⇒ y(t) =− 1 t+c AW: y0=−1

c >0

=⇒ c=−1 y0

<0

=⇒ y(t) = 1

1 y0−t

t y

t^∗ y(t) = 1¹

y0−t

– 81–

(II) Asymptotisch periodische L¨osung

Falls die Zustandsdimensionn= 2 ist muss, ansonsten kanny(t) sich asymptotisch einer periodischen Lösungy_∗(t) nähern, für die gilt

y_∗(t+T) =y_∗(t) f¨ur allet>0 und feste PeriodeT.

Beispiel: siehe obigesLineares Beispiel f¨ur Euler

(III) Chaotisches Verhalten

Falls Dimensionn>2 (einschliesslichn= 2 im nichtautonomen Fall) kann die Lösungy(t) der ODE sich chaotisch verhalten, d.h. auch nach sehr langer Zeit lässt sich keine periodische oder stationäre Struktur erkennen.

Beispiel: Lorenz - Attraktor ( ¨Ubung 2)

– 82–

(2)

Mathematik f¨ur Informatiker III Interpolation mit Polynomen und Splines

Interpolation mit Polynomen (Whd. 1.Semester)

D - 7 Interpolation mit Polynomen und Splines Interpolation mit Polynomen (Whd. 1.Semester)

Satz D.27 (Lagrange - Interpolation)

SeiR=Roder ein anderer K¨orper. Dann gilt:

(i)Es existiert zu jeder Familie von Wertepaaren(xi,yi)∈ R × Rf¨ur i= 0,1, . . . ,n mit unterschiedlichen “Abzissenwerten” xi6=xjf¨ur i6=j ein Interpolationspolynom P(x)vom Grad≤n, so daß

P(xi) =yi f¨ur i= 0,1, . . . ,n.

(ii)Dieses Polynom ist eindeutig und l¨aßt sich darstellen als P(x) =

Xn

i=0

yi (x−x0). . .(x−xi−1)(x−xi+1). . .(x−xn) (xi−x0). . .(xi−xi−1)(xi−xi+1). . .(xi−xn)

| {z }

≡Pi(x)

(iii)Insbesondere folgt aus yi= 0f¨ur i= 0, . . . ,n, dass alle

Koeffizienten ciin P(x) =c0+c1x+c2x²+. . .verschwinden, d.h.

es gilt ci= 0 f¨ur i= 0, . . . ,n.

– 83–

Beispiel – Lagrangepolynom

xi 0 1 2 3

yi -1 2 1 0

P(x) = −1·(x−1)(x−2)(x−3) (0−1)(0−2)(0−3) + 2·(x−0)(x−2)(x−3) (1−0)(1−2)(1−3) + 1·(x−0)(x−1)(x−3) (2−0)(2−1)(2−3) + 0·(x−0)(x−1)(x−2) (3−0)(3−1)(3−2)

P(x) = 2

3x³−4 x² + 19 3 x−1

1 2 3

1 2

−1

– 84–

Warnung:

Interpolationspolynome h¨oherer Ordnung k¨onnen zwischen den vorgegebenen Datenpunktensehr stark oszillieren, deshalb wendet man in der Numerik lieber aus Polynomen niederer Ordnung zusammengesetzte Funktionsmodelle an. =⇒Cubic Splines, Finite Elemente.

8

4

-4 6

2

x

8 6

2 0

10

-2 0

4 PSfrag replacements

Lagrange - Polynom

Kubischer Spline

– 85–

Interpolation durch kubische Splines

Gegeben:

gemessene Datenpaare (xi,yi),i= 0, . . . ,n.

Gesucht:

manipulierbare FunktionP(x) mitP(xi) =yi,i= 0, . . . ,n.

Ansatz

Definiere die interpolierende FunktionP: [x0,xn]→Rin jedem Teilintervall [xi−1,xi] als kubisches PolynomPi, so dass f¨urxi−1≤x≤xi

gilt:

P(x) =Pi(x) =ai(x−xi−1)³+bi(x−xi−1)²+ci(x−xi−1) +di, wobei die 4nKoeffizienten (ai,bi,ci,di) f¨uri= 1, . . . ,nzu bestimmen sind.

– 86–

(3)

Eigenschaften kubischer Polynome

Pihat 4 freie Parameter und die Ableitungen

Pi⁰(x) = 3ai(x−xi−1)²+ 2bi(x−xi−1) +ci

Pi⁰⁰(x) = 6ai(x−xi−1) + 2bi

Pi⁰⁰⁰(x) = 6ai

P⁰⁰⁰⁰i (x) = 0

F¨ur die Bestimmung der 4nKoeffizienten (ai,bi,ci,di),i= 1, . . . ,n, des gesuchten kubischen SplinesP(x) sind genauso viele Gleichungen n¨otig.

Diese werden aus vier verschiedenen Bedingungen, die die interpolierenden Polynome erf¨ullen m¨ussen, hergeleitet.

– 87–

Interpolationsbedingung

Pi(xi) = Pi+1(xi) =yi, i= 1, . . . ,n−1 P1(x0) = y0

Pn(xn) = yn

Mit ∆xi=xi−xi−1folgt aus der Interpolationsbedingung f¨uri= 1, . . . ,n di=yi−1=Pi(xi−1)

ai∆xi³+bi∆xi²+ci∆xi+di=yi=Pi(xi) . Das sindnmal 2 lineare Gleichungen in jeweils 4 Unbekannten.

– 88–

Steigungsbedingung

Pi⁰(xi) =Pi+1⁰ (xi), i= 1, . . . ,n−1 Daraus folgen dien−1 weiteren Bedingungen:

3ai∆xi²+ 2bi∆xi+ci=ci+1, i= 1, . . . ,n−1

Es bleiben nochn+ 1 Freiheitsgrade nach Erf¨ullung der bisher gefundenen 3n−1 linearen Gleichungen.

Kr¨ummungsbedingung

P⁰⁰i(x) =Pi+1⁰⁰ (x), i= 1, . . . ,n Daraus folgenn−1 weitere Bedingungen der Form

6ai∆xi+ 2bi= 2bi+1, i= 1, . . . ,n.

– 89–

Insgesamt hat man nun 4n−2 lineare Gleichungen in 4nUnbekannten, die fehlenden 2 Gleichungen werden durch spezielle Forderungen anP1

undPnim Anfangspunktx0bzw. Endpunktxnerhalten. Diese beiden Bedingungen unterscheiden auch verschiedene Typen kubischer Splines:

Nat¨urlicher kubischer Spline

P⁰⁰(x0) =P⁰⁰1(x0) = 0 =P⁰⁰n(xn) =P⁰⁰(xn)

Im Falle nat¨urlicher Splines sind die letzten fehlenden Gleichungen also b0= 0 und 3an∆xn+bn= 0

Periodischer kubischer Spline

P1(x0) =Pn(xn), P₁⁰(x0) =Pn⁰(xn), P₁⁰⁰(x0) =Pn⁰⁰(xn).

– 90–

(4)

Berechnung der Koeffizienten bei nat¨urlichen Splines

Gesamtbilanz

Man erh¨alt ein sehr strukturiertes lineares Gleichungssystem von 4n Gleichungen in ebenso vielen Unbekannten.

Reduktion auf ein lineares System in (

n−

1) Variablen

zi = P_i+1⁰⁰ (xi) = 2bi+1 f¨uri= 1, . . . ,n−1 z0 = P₁⁰⁰(x0) = 0

zn = Pn⁰⁰(xn) = 0

Lemma D.28

Aus(yi−1,yi,zi−1,zi)ergeben sich die Koeffizienten(ai,bi,ci,di)von Pi

als

di = yi−1 bi = zi−1/2

ai = ^zⁱ_6∆x⁻^zⁱ⁻¹_i ci = ^yⁱ⁻_∆x^yⁱ⁻¹_i −¹6(zi+ 2zi−1) ∆xi – 91–

Struktur des reduzierten Systems bei nat¨urlichen Splines

Mit

αi= 2(∆xi+ ∆xi+1) und βi= ∆xi

sowie

ri= 6 yi+1−yi

∆xi+1 −yi−yi−1

∆xi

ist zur Bestimmung derzi,i= 1, . . . ,n−1, das folgende

diagonaldominante symmetrische tridiagonale lineare Gleichungssystem zu l¨osen:







α1 β2

β2 α2 β3

β3 α3 β4

. .. ... . ..

βn−2 αn−2 βn−1

βn−1 αn−1











 z1

z2

z3

... zn−2

zn−1







=





 r1

r2

r3

... rn−2

rn−1







– 92–

Mathematik f¨ur Informatiker III Numerische Integration – Quadratur

D - 8 Numerische Integration – Quadratur

Gr¨unde f¨ur numerische Integration

IFunktionen ohne geschlossen darstellbare Stammfunktion

IStammfunktion nur durch sehr komplizierte Formel darstellbar Beispiele D.29 (Funktionen ohne geschlossenes Integral)

IR

e⁻^x²dx Gauß’sche Glockenkurve

IR√

1−k²sin²t dt Elliptisches Integral

– 93–

Interpolatorische Quadraturformeln

Um eine N¨aherung des bestimmten Integrals Zb

a

f(x)dx

zu berechnen, wird das Integrationsintervall [a,b] inn∈IIgleichgrosse Teilintervalle [x0,x1], . . . ,[xn−1,xn] der Längehn=^b⁻_nâunterteilt. Dabei giltxi=a+i∗hnund insbesonderex0=aundxn=b. Mitfi=f(xi) wird der Funktionswert an deri-ten Stützstelle bezeichnet.

Riemann’sche Summen

Zb

a

f(x)≈ Xn

i=1

f(xi)hn= Xn

i=1

fihn

Fehlerterm Riemann’sche Summen

Zb a

f(x)− Xn

i=1

fihn

≤b−a

2 ·hn · max

x∈[a,b]|f⁰(x)|

– 94–

(5)

Summierte Trapezregel

Tn=hn

"

1 2(f0+fn) +

n−1

X

i=1

fi

#

Approximationsfehler summierte Trapezregel

Zb

a

f(x)dx−Tn

≤b−a

12 ·h²n · max

x∈[a,b]|f⁰⁰(x)|

– 95–

Kepler’sche Fassregel

Ansatz: Quadratischer Splineg(x) durch die Punkte (a,f(a)),(^a+b₂,f(^a+b₂)), und (b,f(b))

g(x) =cx²+dx+e Durch geeignete Umformung des Ansatzes erh¨alt man eine

Berechnungsvorschriftohnedie Koeffizientenc,dundedes Splinesg(x):

S0=b−a 6

f(a) + 4f(^a+b₂ ) +f(b)

– 96–

Simpson’sche Regel (Summierte Kepler’sche Fassregel)

Anwendung der Fassregel auf die Teilintervalle der L¨angehn=^b⁻n^a, n gerade, ergibt die Simpson’sche Regel:

Sn=hn

3



f0+fn+ 2

n 2−1

X

i=1

f2i+ 4 Xn/2

i=1

f2i−1





Approximationsfehler summierte Simpson’sche Regel

Zb

a

f(x)dx−Sn

≤b−a

180 ·hn⁴· max

x∈[a,b]|f⁽⁴⁾(x)|

– 97–

Quadratur mit Extrapolation – Romberg’s Verfahren

F¨ur hinreichend oft differenzierbare Integrandenf(x) beschreibt die Euler-Maclaurinsche Summenformelden Fehler der summierten TrapezregelTnals Polynom in geraden Potenzen der Schrittweitehn:

Tn= Zb

a

f(x)dx+ XN

k=1

α2kh^2kn +O(h^2N+2n )

Die dabei auftretenden Koeffizientenα2ksind vonhnunabh¨angige Konstanten.

Damit k¨onnen Fehlerterme von Quadraturformeln durch sog.

Extrapolation zur Grenze/zum Limiteliminiert werden, in der Werte einer Quadraturformel bei unterschiedlichen Schrittweitenhn,n=n1,n2, . . ., kombiniert werden.

Bei geschickter Wahl der Extrapolation erreicht man eine Aufhebung von Fehlertermen kleiner Ordnung, so das der extrapolierte Wert eine deutlich genauere Approximation des gesuchten Integralwertes ist.

– 98–

(6)

Romberg Verfahren

Zuerst wird f¨urn= 1 die Trapezregel auf dem gesamten Integrationsintervall [a,b] ausgewertet. Der erhaltene WertT1(d.h.

Schrittweiteh1=b−a) wird als erster EintragR₀⁰in die erste Zeile der Tabelle eingetragen.

Mit halbierter Schrittweiteh2=h1/2 wirdT2=R1⁰berechnet und in die erste Spalte der zweiten Zeile direkt unterR0⁰notiert:

k n= 2^k R_k⁰

0 1 R0⁰

1 2 R₁⁰ R₁¹ Daraus berechnet man denextrapolierten WertR1¹mittels

R₁¹=4R₁⁰−R₀⁰

3 = S2,

was aber genauSimpsons Regelf¨urn= 2 ergibt.

– 99–

Romberg Verfahren (Fortsetzung)

Dieses Vorgehen kann in einer neuen Zeile der Tabelle fortgef¨uhrt werden.

Diek-te Zeile erhält man dabei, indem zunächst die Trapezregel mit erneut halbierter Schrittweitehn=h2^k(d.h.n= 2^k) ausgeführt wird und T2^kalsR⁰kin die erste Spalte eingetragen wird.

In den darauffolgendenkExtrapolationsschritten werden jeweils die WerteR_k^jderk-ten Zeile f¨urj= 1, . . . ,kaus dem links stehenden Wert Rk^j⁻¹und dem links dar¨uber stehenden WertRk^j⁻−1¹berechnet:

R_k^j=4^jRk^j⁻¹−Rk^j⁻−1¹

4^j−1 =R_k^j⁻¹+ 1 4^j−1

R_k^j⁻¹−R_k^j⁻¹₋₁

j= 1, . . . ,k Insgesamt ergibt sich damit das folgende Tableau:

– 100–

Romberg Verfahren (Fortsetzung)

k n= 2^k Rk⁰=Tn Rk¹ R²k . . . 0 1 R₀⁰=T1

1 2 R₁⁰=T2 R₁¹ 2 4 R2⁰=T4 R2¹ R²2

3 8 R₃⁰=T8 R₃¹ R²₃ R₃³ 4 16 R⁰₄=T16 R₄¹ R²₄ R₄³ R₄⁴

... ... ... ... . ..

AlsAbbruchbedingungeignet sich die Differenz zwischen den beiden zuletzt berechneten Diagonalelementen des Schemas.

Falls mit einer vorgegebenen Gr¨osseδdie Bedingung

|R^kk−R^kk−1⁻¹| ≤δ

erf¨ullt ist, dann wird das Verfahren beendet undRk^kals N¨aherung des IntegralsRb

af(x)dxbetrachtet. – 101–

Approximationsfehler Romberg-Verfahren

F¨urf∈C^2k+2([a,b]) gilt:

Rk^k− Zb

a

f(x)dx

≤(b−a)h²₁h²₂. . .h²₂kα2k+2 max

x∈[a,b]|f^(2k+2)| wobeiα2k+2wiederum eine Konstante ist.

Bemerkung

Die auftretenden Konstantenαiergeben sich als αi=Bi

i!,

wobei dieBidie so genanntenBernoulli - Zahlensind. Diese berechnen sich rekursiv aus

Bi= (−1)ⁱ⁻¹

"

2i−1 2(2i+ 1)+ (2i)!

i−1

X

k=1

Bk

(2i−2k+ 1)!(2k)!

# .

– 102–