HUMBOLDT–UNIVERSIT ¨ AT ZU BERLIN

(1)

HUMBOLDT–UNIVERSIT ¨ AT ZU BERLIN

Mathematisch-Naturwissenschaftliche Fakult¨at II Institut f¨ur Mathematik

Prof. PhD. Andreas Griewank Jan Riehme

Humboldt-Universit¨at zu Berlin, Institut f¨ur Mathematik, Unter den Linden 6, D-10099 Berlin

Musterl¨ osung ¨ Ubungsserie 2 Mathematik f¨ ur Informatiker III

Aufgabe 1: L¨ose folgende Anfangswertaufgaben:

(i) y⁰= (y³+1)sin(x)/y², y(0) =1 L¨osung:Die ODE

y⁰ = (y³+ 1)sin(x)/y²=y³+ 1 y²

| {z }

g(y)

sin(x)

| {z }

f(x)

ist eine separable Differentialgleichung erster Ordnung.

Verwende L¨osungsmethode auf Folie nach Satz C.12:

F(x) :=

Z x x₀

f(t)dt= Z x

0

sin(t)dt= −cos(t)

x

0 = 1−cos(x) und

G(y) :=

Z y y₀

1 g(t) dt=

Z x 0

t²

t³+ 1 dt= 1

3 ln|t³+ 1|

y 1 =1

3 ln|y³+ 1| −ln 2

=:z(y) Nun UmkehrfunktionG⁻¹(z) =y(z) bilden:

ln|y³+ 1|= 3z+ ln 2 =⇒ |y³+ 1|=e^{3z+ln 2}=e^3ze^{ln 2}= 2e^3z Fallunterscheidung wegen Betragsfunktion:

Fall 1: y³≥ −1 =⇒ |y³+ 1|=y³+ 1≥0

=⇒y³= 2e^3z−1 =⇒y=√³

2e^3z−1

Einsetzen vonF(x) als Argument vonG⁻¹(z) liefert L¨osung der ODE y1(x) =p³

2e^{3−3 cos(x)}−1 Pr¨ufung der Anfangsbedingung

y₁(0) =p³

2e^{3−3 cos(0)}−1 = p³

2e³⁻³−1 = √³

2−1 =√³

1 = 1 =y₀ Damit isty1(x)eine L¨osungdes gegebenen Anfangswertproblems.

Fall 2: y³<−1 =⇒ |y³+ 1|=−(y³+ 1)>0

=⇒y³=−2e^3z−1 =⇒y=√³

−2e^3z−1 =−√³

2e^3z+ 1

Einsetzen vonF(x) als Argument vonG⁻¹(z) liefert ( f¨ury³<−1 ) y2(x) =−p³

2e^{3−3 cos(x)}+ 1 Pr¨ufung der Anfangsbedingung

y₂(0) =−p³

2e^{3−3 cos(0)}+ 1 =−p³

2e³⁻³+ 1 =−√³

2 + 1 =−√³

36= 1 =y₀ Damit isty2(x)keine L¨osungdes gegebenen Anfangswertproblems.

(2)

(ii) 2(1+x) +3t ˙x=0, x(1) =0 oder x(1) =−2 L¨osung:Hier gibt es (mindestens) zwei L¨osungsvarianten:

Variante A: L¨osung als Separable ODE

˙ x= −1

3t

|{z}

f(t)

(1 +x

| {z }

g(x)

)

Verwende L¨osungsmethode auf Folie nach Satz C.12:

F(t) :=

Z t t₀

f(y)dy= Z t

1

−2

3y dy=−2 3

Z t 1

1

y dy=−2 3 ln|y|

t 1

=−2 3ln|t|

und, wenn wir t >= 1 =t₀ vorausetzen (dieser Fall reicht uns hier), F(t) =−2

3lnt,

sowie (Achtung: Hier wird der Anfangswertx0 nicht sofort eingesetzt, da in dieser Musterl¨osung alle Anfangswerte diskutiert werden sollen. Wenn nur ein nur ein Anfangswert vorliegt, dann sollte nat¨urlich dieser gleich eingesetzt werden.)

G(x) :=

Z x x₀

1 g(y) dy=

Z x x₀

1

1 +y dy= ln|y+ 1|

x

x₀= ln|x+ 1| −ln|x0+ 1|=:z(x) Nun UmkehrfunktionG⁻¹(z) =x=x(z) bilden:

ln|x+ 1|=z+ ln|x0+ 1| =⇒ e^z+ln^|x⁰^+1|= e^z· |x0+ 1|=|x+ 1|

Fallunterscheidung wegen Betragsfunktion f¨urx:

Fall 1: x≥ −1 =⇒ |x+ 1|=x+ 1≥0 =⇒x= x(z) =e^z|x0+ 1| −1 Einsetzen vonF(t) =−²₃lntals Argument von G⁻¹(z) liefert (x≥ −1 )

x=e⁻²³^ln^t|x0+ 1| −1 = e^ln^t⁻²3|x0+ 1| −1 = |x0+ 1|

√3

t² −1 =x1(t) Fall 2: x≤ −1 =⇒ |x+ 1|=−(x+ 1)≥0 =⇒x= x(z) =−e^z|x0+ 1| −1

Einsetzen vonF(t) =−²₃lntals Argument von G⁻¹(z) liefert (x <−1 )

−e⁻²³^ln^t|x0+ 1| −1 =− e^ln^t⁻²₃

|x0+ 1| −1 = −|x0+ 1|

√3

t² −1 =x₂(t)

Einsetzen der verschiedenen Anfangsbedingungen (nurx0wird jetzt verwendet,t0wurde schon bei F(t)benutzt):

x0=0

x₁(t) = ^|x^√⁰₃^+1|

t² −1 = ^|0+1|^√₃

t² −1 = ^√₃¹

t² −1 = t^−2/3−1=x₁(t) Pr¨ufe Anfangsbedingung:x1(1) = 1^−2/3−1 = 1−1 = 0 =x0

Also istx1(1) =t^−2/3−1eine L¨osung des Anfangswertproblems mitx(1) = 0.

x2(t) =−^|x√3⁰^+1|

t² −1 =−^|0+1|√3

t² −1 =−√3¹

t² −1 = −t^−2/3−1=x2(t) Pr¨ufe Anfangsbedingung:x2(1) =−1^−2/3−1 =−1−1 =−26= 0 =x0

Also istx2(1) =−t^−2/3−1 keine L¨osungdes Anfangswertproblems mitx(1) = 0.

x0=−2

x1(t) = ^|x√⁰3^+1|

t² −1 = ^|−2+1|√3

t² −1 = √3¹

t² −1 = t^−2/3−1=x1(t) Pr¨ufe Anfangsbedingung:x1(1) = 1^2/3−1 = 1−1 = 06=−2 =x0

Also istx1(1) =t^−2/3−1keine L¨osungdes Anfangswertproblems mitx(1) =−2.

x2(t) =−^|x√3⁰^+1|

t² −1 =−^|−2+1|√3

t² −1 =−√3¹

t²−1 = −t^−2/3−1=x2(t) Pr¨ufe Anfangsbedingung:x2(1) =−1^−2/3−1 =−1−1 =−2 =x0

Also istx2(1) =−t^−2/3−1 eine L¨osungdes Anfangswertproblems mitx(1) =−2.

(3)

x0=−1 (das ist dasx0 der urspr¨unglichen Aufgabe) x1(t) = ^|x√⁰3^+1|

t² −1 = ^|−1+1|√3

t² −1 = 0−1 = −1=x_∞=x1(t) x₂(t) =−^|x√3⁰^+1|

t² −1 =−^|−1+1|√3

t² −1 =−0−1 = −1=x_∞=x₂(t)

Der Anfangswertx₀=−1 führt auf die Konstante−1 als Lösung, das ist aber genau der vom Anfangswert unabhängige asymptotische Grenzwertx_∞ der Lösung der ODE!

Dieses Verhalten kann wie folgt interpretiert werden: Der zugrunde Prozess (z.Bsp. eine che- mische Reaktion) strebt immer einem Gleichgewichtszustand entgegen, den es u.U. nur im Grenzwert für unendlich lange Prozesszeiten erreicht. Startet man den Prozess allerdings schon im Gleichgewicht ( alsox?−=X∞), dann bleibt das Gleichgewicht erhalten. (Zumin- dest in so einfachen Modellen, bei denen keine anderen äußeren Einflüsse eine Rolle spielen).

0 1 2 3 4 5 6 7 8 9 10 t

x(t)

−2

−3

x_∞=−1

x(t) =t^−2/3−1

x(t) =−t^−2/3−1

Variante B: L¨osung als lineare ODE erster Ordnung

2(1 +x) + 3tx˙ = 0, x(1) = 0 oderx(1) =−2 Umschreiben:

˙ x+ 2

3t

|{z}

a(t)

x=−2 3t

|{z}

f(t)

Schritt 1: Finde L¨osungx_h(t) der homogenen ODE (SatzD.15)

˙ x+ 2

3tx= 0 ^Ansatz=⇒ x(t) =ce^−A(t) mit A(t) eine Stammfunktion vona(t)

A(t) = Z

a(t)dt= 2 3

Z 1 t dt=2

3ln|t|

Also

xh(t) =c e^−A(t)=c e⁻²³^ln^|t|=c

e^ln^|t|−²₃

=c t⁻²³ = c

√3

t² =xh(t) Schritt 2: Finde spezielle L¨osung der inhomogenen ODE mit dem Ansatz aus Satz D.16:

x(t) = Z t

t0

f(z)e^A(z)dz·e^−A(t)+xh(t) = Z t

t0

f(z)e^A(z)dz+c

e^−A(t)

= Z t

t₀

−2

3ze²³^ln^|z|dz+c

t⁻²³ =

−2 3

Z t t₀

1

zz²³dz+c

t⁻²³

=

−2 3

Z t t₀

z⁻¹³dz+c

t⁻²³ = −2 3

"

3 2z²³

t

t0

# +c

! t⁻²³

=

−62 63·63

62 h

t²³ −t

2 3

0

i +c

t⁻²³ =

t₀=1

−t²³ + 1 +c t⁻²³

= (c+ 1)t⁻²³ −t²³ ·t⁻²³ = (c+ 1)t⁻²³ −1 =x(t)

Einsetzen der verschiedenen Anfangsbedingungen, um die verbliebene Konstante c zu bestimmen (nurx₀ wird jetzt verwendet, t₀ wurde schon beix(t)benutzt):

(4)

x(1) =0

x(1) = 0 = (c+ 1)1⁻²³ −1 =c+ 1−1 =c= 0 =⇒ c= 0 Also ist

x(t) =t⁻²³ −1 = 1

√3

t² −1 eine L¨osung der ODE mit der Anfangsbedingung x(1) = 0.

Pr¨ufe Anfangsbedingung:x(1) = 1^−2/3−1 = 1−1 = 0 x(1) =−2

x(1) =−2 =c =⇒ c=−2 Also ist

x(t) =−t⁻²³ −1 =− 1

√3

t² −1 eine L¨osung der ODE mit der Anfangsbedingung x(1) =−2.

Pr¨ufe Anfangsbedingung:x(1) =−1^−2/3−1 =−1−1 =−2 x(1) =−1 (das ist dasx₀der urspr¨unglichen Aufgabe)

x(1) =−1 =c =⇒ c=−1 Also ist

x(t) = (−1 + 1)t⁻²³ −1 = 0·t⁻²³ −1 =−1 eine L¨osung der ODE mit der Anfangsbedingung x(1) =−1.

Der Anfangswertx0=−1 führt auf die Konstante−1 als Lösung, das ist aber genau der vom Anfangswert unabhängige asymptotische Grenzwertx_∞ der Lösung der ODE!

(iii) y⁰⁰+y⁰+y=0, y(0) =1=y⁰(0) L¨osung:

Das ist eine lineare ODE zweiter Ordnung mit konstanten Koeffizienten (siehe Definition D.17):

Die Koeffizienten sinda1=a2= 1.

Schritt 1:Suche L¨osung der homogenen ODE, dazu stellen wir das charakteristische Polynom auf und betrachten dessen Wurzeln:

0 =λ²+λ+ 1 λ_1,2=−1

2± r1

4 −1 =−1 2 ±

√3 2 i Somit ist λ = −¹₂ +

√3

2 i (und ihre konjugiert komplexe Zahl ¯λ = −¹₂ −

√3

2 i ) eine einfache komplexe Nullstelle des charakteristischen Polynoms. Also sind

e⁻¹²^xcos^√

3 2 x

und e⁻¹²^xsin^√

3 2 x

L¨osungen der homogenen linearen ODE ( siehe Folie 60 nach Definition D.20 ).

Schritt 2: Konstruktion der L¨osung der Ausgangs - ODE (nicht homogen)

Nach Satz D.19 spannen dien= 2 Lösungen der homogenen linearen ODEn-ter Ordnung einen Vektorraum der Dimensionn= 2 auf. Die Lösungen der homogenen ODE bilden also eine Basis des Vektorraumes der inhomogenen Lösungen.

Eine L¨osung der inhomogenen ODE gewinnen wir aus dem Ansatz ( Linearkombination der L¨osungen der homogenen ODE mit den Koeffizientenc1 undc2)

y(x) =c1·e⁻¹²^xcos^√

3 2 x

+c2·e⁻¹²^xsin^√

3 2 x

=e⁻¹²^xh

c1·cos^√

3 2 x

+c2·sin^√

3 2 xi

und bestimmen mit Hilfe der beiden Anfangsbedingungen y(0) = 1 = y⁰(0) die Konstanten c₁ undc₂:

y(0) = 1 =e⁻¹²^·0h

c1·cos^√

3 2 ·0

+c2·sin^√

3 2 ·0i

= 1 [c1·1 +c2·0] = c1= 1

(5)

Zur Bestimmung des letzten Koeffizienten verwenden wir y⁰(0) = 1, m¨ussen also erst einmal die erste Ableitung des L¨osungansatzesy(x) =e⁻¹²^xh

c₁·cos^√

3 2 x

+c₂·sin^√

3 2 xi

bestimmen:

y⁰(x) = e⁻¹²^xh

c1·cos^√

3 2 x

+c2·sin^√

3 2 xi0

=−1 2e⁻¹²^xh

c1·cos^√

3 2 x

+c2·sin^√

3 2 xi + e⁻¹²^xh

−c1·sin^√

3 2 x

·

√ 3

2 +c2·cos^√

3 2 x

·

√ 3 2

i

Mity⁰(0) = 1 undc1= 1 erhalten wir y⁰(0) = 1 =−1

2e⁰[c1·cos 0 +c2·sin 0] + e⁰h

−c1·

√3

2 ·sin 0 +c2·

√3 2 ·cos 0i

=−1

2 c1 + h c2·

√3 2

i

c₁==1c2·

√3 2 −1

2

= 1!

und somit

c2= 3 2· 2

√3 =√ 3.

Insgesamt ergibt sich damit als L¨osung der inhomogenen ODE 2.ter Ordnung y(x) =e^−x/2 h√

3 sin^√

3 2 x

+ cos^√

3 2 xi

Aufgabe 2: Numerische L¨osung von Anfangswertproblemen

Wende den expliziten Euler und die Mittelpunktsregel auf das Problem

˙

y=y+t, y(0) = 1, T = 1

an. Verwende dabei die Schrittweiten h= 1/nk mitnk = 2^k f¨urk= 1,2,3, . . . ,9.

Berechne die exakte L¨osung y(1) und plotte f¨ur beide Methoden −log|ynk−y(1)| uber¨ k in einem Diagramm. Hierbei ist y_n_k der von der Methode in n= 2^k Schritten erhaltene Wert.

L¨osung:

Bestimme exakte L¨osung der linearen ODE 1.Ordnung (Def. D.14):

˙

y−y = ˙y+ (−1

|{z}

a(t)

)y= t

|{z}

f(t)

• L¨ose homogene ODE ˙y−y= 0 (Satz D.15):

– BestimmeA(t) als Stammfunktion vona(t) =−1 =⇒A(t) =−t – L¨osung der homogenen ODE ist damityh(t) =c·e^−A(t)=c·e^t

• L¨osung der inhomogenen ODE (Satz D.16) y(t) =

Z t t0

f(z)e^A(z)dz+c

·e^−A(t)= Z t

0

ze^−zdz+c

·e^t

Partielle Integration [R

uv⁰dx=uv−R

u⁰vdx] mitu=z,v⁰=e^−z,v=−e^−z ergibt:

Z t 0

z e^−zdz=−z e^−z

t

0

− Z t

0

−e^−zdz= −z e^−z−e^−z

t

0

= (−t−1)e^−t−(−1) = 1−(t+ 1) e^t und damit

y(t) =

1−t+ 1 e^t +c

e^t=e^t−(t+ 1) +c e^t= (1 +t)e^t−t−1

• Anfangswerty(0) = 1 benutzen, um Konstatneczu bestimmen y(0) = 1 =

1−e⁻⁰(0 + 1)

e⁰= 1−1 +c= 1 =⇒ c= 1

(6)

Insgesamt

y(t) = 2e^t−t−1, also y(1) = 3.4365636569180904707

Ergebnisse von Eulerverfahren und Mittelpunkt-Regel im Vergleich mit exakter L¨osung:

k y₂k,Eul |y₂k,Eul−y(1)| y₂k,M id |y₂k,M id−y(1)|

1 2.500000 0.936564 3.281250 0.155314

2 2.882812 0.553752 3.389711 0.046853

3 3.131569 0.304995 3.423682 0.012882

4 3.275857 0.160707 3.433187 0.003377

5 3.353980 0.082584 3.435699 0.000865

6 3.394690 0.041874 3.436345 0.000219

7 3.415478 0.021086 3.436509 5.5e-05

8 3.425983 0.010581 3.436550 1.4e-05

9 3.431264 0.0053 3.436560 4e-06

Tabelle −log|ynk−y(1)|:

k −log|yn_k−y(1)|

Euler Midpoint 1 0.065538 1.862309 2 0.591040 3.060756 3 1.187461 4.351988 4 1.828175 5.690882 5 2.493947 7.053580 6 3.173096 8.428112 7 3.859164 9.808536 8 4.548751 11.191898 9 5.240113 12.576727

-2 0 2 4 6 8 10 12 14

0 2 4 6 8 10 12

Euler Midpoint

Abbildung 1:−log|yn_k−y(1)|f¨ur Euler und Mittelpunkt

M¨oglichesC- Programm:

1 #include <math . h>

2

3 #define MAXK 1 0

(7)

4

5 // r e c h t e S e i t e d e r ODE

6 double f (double t , double y ) { 7 return y + t ;

8 } 9

10 double T = 1 . 0 ; 11 double T0 = 0 . 0 ; 12 double Y0 = 1 . 0 ; 13

14 i n t main ( void ) { 15 long i n t k , n , i ;

16 double h , tk , yk , yk0 , f k ; 17 double ykm ;

18

19 n=1;

20 f o r ( k = 1 ; k <MAXK; k ++ ) {

21 n ∗= 2 ;

22 h = (T−T0 ) / (double) n ;

23 p r i n t f ( ”#\n# k=%d n=%d h=%l f\n#\n” , k , n , h ) ; 24

25 t k = T0 ;

26 yk = ykm = Y0 ;

27 p r i n t f ( ” % l f % l f % l f % l f\n” ,

28 tk , yk , ykm , 2 . 0∗exp ( t k )−tk−1 ) ; 29 f o r ( i = 1 ; i < n + 1 ; i ++ ) {

30 // E u l e r

31 yk += h ∗ f ( tk , yk ) ;

32 // m i d p o i n t

33 f k = f ( tk , ykm ) ;

34 ykm += h ∗ f ( ( h∗( (double) i−. 5 ) ) , ykm + h∗f k / 2 . 0 ) ;

35 // t i m e s t e p

36 t k = i∗h ;

37 p r i n t f ( ” % l f % l f % l f % l f\n” ,

38 tk , yk , ykm , 2 . 0∗exp ( t k )−tk−1 ) ;

39 }

40 p r i n t f ( ”\n\n\n” ) ;

41 }

42 }

Aufgabe 3: Numerische L¨osung von Systemen von ODEs

Vom einem Startpunkt nahe am Ursprung (z.B.(0.1,0.1,0.1)) wende die Mittelpunktsregel auf das folgende System von drei gew¨ohnlichen Differentialgleichungen erster Ordnung an:

dx

dt = −10x+ 10y dy

dt = 28x−y−xz dz

dt = −8/3z+xy.

L¨osung f¨ur T = 1, k = 1, . . . ,35 mit a_k = (x_k, y_k, z_k)^T mit x_k, y_k und z_k jeweils berechnet mit Schrittweite 2^−k:

k x y z kak−ak−1k

1 140.495916 -43.980152 -31.397437 —-

2 305.485175 213.513767 -225.221352 3.6e+02

3 9481751.390046 35825683761.798683 35446260109.048477 5.0e+10

4 3.990979 0.828050 33.591529 5.0e+10

5 -8.085631 -8.340393 26.680917 1.7e+01

(8)

k x y z kak−a_k−1k 6 -8.662505 -9.682238 26.273605 1.5e+00 7 -8.694260 -9.813123 26.166494 1.7e-01 8 -8.695994 -9.825330 26.152600 1.9e-02 9 -8.695807 -9.825945 26.150978 1.7e-03 10 -8.695686 -9.825797 26.150809 2.6e-04 11 -8.695647 -9.825722 26.150796 8.5e-05 12 -8.695636 -9.825699 26.150796 2.6e-05 13 -8.695633 -9.825692 26.150797 7.1e-06 14 -8.695633 -9.825691 26.150797 1.8e-06 15 -8.695632 -9.825690 26.150797 4.7e-07 16 -8.695632 -9.825690 26.150797 1.2e-07 17 -8.695632 -9.825690 26.150797 3.0e-08 18 -8.695632 -9.825690 26.150797 7.5e-09 19 -8.695632 -9.825690 26.150797 1.9e-09 20 -8.695632 -9.825690 26.150797 4.7e-10 21 -8.695632 -9.825690 26.150797 1.2e-10 22 -8.695632 -9.825690 26.150797 3.3e-11 23 -8.695632 -9.825690 26.150797 5.0e-12 24 -8.695632 -9.825690 26.150797 7.7e-12 25 -8.695632 -9.825690 26.150797 7.4e-12 26 -8.695632 -9.825690 26.150797 4.9e-12 27 -8.695632 -9.825690 26.150797 4.6e-12 28 -8.695632 -9.825690 26.150797 2.4e-11 29 -8.695632 -9.825690 26.150797 1.6e-11 30 -8.695632 -9.825690 26.150797 3.3e-11 31 -8.695632 -9.825690 26.150797 4.6e-11 32 -8.695632 -9.825690 26.150797 1.6e-10 33 -8.695632 -9.825690 26.150797 4.0e-10 34 -8.695632 -9.825690 26.150797 5.0e-10 35 -8.695632 -9.825690 26.150797 6.3e-09

In dieser Tabelle sieht man sehr schnell, das f¨urk <5 durch die zu grobe Diskretisierung mith_k =₂¹_k keine sinnvollen Werte berechnet werde (nicht einmal die Vorzeichen vonxundy sind korrekt). Da ist die Auswirkung des Diskretisierungsfehlers.

Erst abk≥5 ist die Diskretisierung so fein, das der Diskritiesierungsfehler nur noch ein untergeordnete Rolle spielt und die berechneten Werte gegen den zum Zeitpunkt T = 1 gehörigen Punkt auf der Lösungstrajektorie konvergieren. Dabei verringert sich mit jeder Halbierung der Schrittweite auch der Abstand kak−a_k−1k der zugehörigen Lösung a_k = (x_k, y_k, z_k)^T zur vorhergehenden Lösung a_k−1 ( siehe Spalte 5 der obigen Tabelle ).

Abk= 23 fängt der Abstandkak−a_k−1kwieder an zu wachsen. Der Grund dafür sind die nun immer stärker werdenden Rundungsfehler bei der Berechnung vonak mitk≥23.

Verifiziere zunächst über eine kurze Integrationsperiode T = 1 mit den gleichen Zeit- schritten wie in Aufgabe 2, dass die numerische Integration tatsächlich mit der Ordnung zwei erfolgt. Dazu berechne die annähernde Fehlerkonstante

4 3

kan_k−an_k+1k h²_n_k

mit an_k = (xn_k, yn_k, zn_k)^T und untersuche, ob sie im Wesentlichen von k unabh¨angig ist.

Fehlerkonstante ˜c=⁴₃^ka^nk^−a_h2^nk+1^k

nk f¨urk= 1, . . . ,35:

k ˜c k ˜c k ˜c

1 1931.022318 10 119.184565 23 722.843709

2 1075148583592.741089 11 144.195121 24 2784.797539 3 4300594331743.996094 12 157.957615 25 7423.270558

4 5687.697234 13 165.062863 26 27566.738605

5 2070.285777 14 168.663292 27 584342.038577 6 939.793479 15 170.473983 28 1557510.077861 7 405.799112 16 171.381771 29 12709087.405301

(9)

k ˜c k ˜c k c˜ 8 152.447179 17 171.859347 30 70009277.896700 9 89.219573 18 172.113293 31 1000007865.513323 19 171.264984 32 9827056893.592941 20 172.066506 33 49666443572.844818 21 190.881547 34 2470013054559.398438 22 117.834077

Interpretation:

Die Gr¨oßec=c(T), die bei einem Integrator der Ordnungpals Vorfaktor des Diskretisierungsfehlers h^p in

kak−y(T)k = c(T)h^p_k+O(h^p+1) ≈ c(T)h^p_k

auftritt, wird Fehlerkonstante genannt, da sie nicht von der gewählten Schrittweitehk abhängt (aber natürlich von der IntegrationsdauerT).

Eine Näherung vonc(T) für die Mittelpunktsregel erhält man aus Beispiel D.25 des Skripts. Dort wird aus der mit der Mittelpunktsregel und Schrittweitehk berechneten Lösung ak und der mit doppelter Schrittweiteh_k−1berechnetena_k−1 die Näherung ˜c(T) hergeleitet:

c(T)≈⁴₃kak−a_k+1k

h²_k ≡ ˜c(T) = ˜c.

F¨urk= 1, . . . ,4 macht die Betrachtung ˜c keinen Sinn, da zur Berechnung Werte benutzt werden, die mit der L¨osung aufgrund der zu grossen Schrittweite nichts zu tun haben (siehe oben).

Fürk= 5, . . . ,9 ist eine deutliche Reduktion von ˜czu sehen, allerdings sind die Schwankungen immer noch zu stark für eine Aussage über die Abhängigkeit der Fehlerkonstanten von der Schrittweite.

Trotzdem ist bereits hier bemerkenswert, das die ˜c als N¨aherung des Vorfaktors im Fehler c(T)h²_k+ O(h²⁺¹_k ) nicht mit kleiner werdender Schrittweite w¨achst, sondern kleiner wird.

Wird die Schrittweite weiter verkleinert (k = 10, . . . ,22), so kann man an den berechneten Werten für ˜c sehen, dass die Fehlerkonstante tatsächlich von der Wahl der Schrittweite unabhängig ist: Die entsprechenden Werte in der Tabelle 3 schwanken nur noch sehr wenig, auch wenn die Schrittweite sehr stark verkleinert wird ( und zwar um den Faktor 2²²/2¹⁰= 2¹²= 4096).

Sehr sch¨on kann man hier wieder sehen, dass die Mittelpunktsregel ein Verfahren der Ordnungp= 2 ist: Das ˜c = ⁴₃^ka^k^−a_h2^k+1^k

k

trotz des quadratisch wachsenden Wertes von _h¹

k nahezu konstant bleibt impliziert sofort, das die Differenzkak−a_k+1k auch mit Ordnung 2 kleiner wird.

Bei weiterer Verkleinerung der Schrittweite (hier abk= 23, . . . ,35 durchgeführt) nimmt der berechnete Wert der Fehlerkonstanten wieder zu, was aber nicht als Abhängigkeit von der Schrittweite interpretiert werden darf: Die weitere Verkleinerung der Schrittweite führt nun aufgrund immer stärker ins Gewicht fallender Rundungsfehler zu wieder mit k wachsenden Veränderungen kak −ak+1k der Lösung ak

(Spalte 5 der Tabelle 3), die noch durch den Faktor _h¹2

k weiter verst¨arkt werden.

Insgesamt hat man damit experimentell nachgewiesen, das die Fehlerkonstantec(T) nicht vonk(und damit der Schrittweiteh_k) abh¨angig ist.

Dann integriere über den längeren Zeitraum T = 20 (besser T = 25 oder T = 30) und beobachte die Abhängigkeit der Lösung von der Schrittzahl und dem Anfangspunkt, der immer noch nahe am Ursprung variiert werden kann. Stelle dabei sicher, das trotz nun deutlich vergrössertem Zeitintervall die gleichen Zeitschritteh_k = 1/2^k verwendet werden!

Im folgenden wird jeweils die x-Komponente des Anfangswertes a⁰= (x0, y0, z0) = (0.1,0.1,0.1) mit einen ∆ gestört. ak = (xk(T), yk(T), zk(T))^T = (xk, yk, zk)^T bezeichnet dabei die mit Schrittweite hk = 1/2^k berechnete numerische Lösung vom ungestörten Anfangswerta⁰ = (.1, .1, .1), während die zum gestörten Anfangswert ã⁰ = (.1 + ∆, .1, .1) gehörige Lösung mit ãk = (˜xk,y˜k,z˜k)^T bezeichnet wird.

∆ = 10⁻⁴

In der folgenden Tabelle sind die Differenzen zwischenak und ˜ak angeben:

T = 20 T = 25 T = 30

k kak−ãkk kak−ãkk kak−ãkk

1 NaN NaN NaN

2 NaN NaN NaN

3 NaN NaN NaN

(10)

T = 20 T = 25 T = 30 k ka_k−a˜_kk ka_k−a˜_kk ka_k−a˜_kk 4 20.501975 13.305980 7.298038 5 0.002802 0.008084 2.195777 6 0.012778 0.395098 3.518372 7 0.017659 2.864944 26.290506 8 0.009686 28.938541 19.557223 9 0.012112 11.627949 15.387764 10 0.014431 0.781808 12.846552 11 0.015420 1.312642 3.071959 12 0.015722 3.318790 13.308455 13 0.015804 9.672175 12.328687 14 0.015825 13.202092 16.299790 15 0.015830 14.267615 18.440682

Die Schrittweiten hk = 1/2^k für k = 1,2,3 sind für die Integration über den langen Zeitraum von T = 20 (und größere Integrationsintervalle) viel zu grob, der Diskretisierungsfehler führt während der Berechnung zu im Gleitkommaraster nicht darstellbaren Zahlen. Im weiteren werden wir deshalb immer mitk= 5 starten.

Fürk >3 beobachtet man, dass fürT = 20 die Auswirkungkak−ãkkder Störung des Startwertes auf die Lösung scheinbar beschränkt ist, während für größere Integrationsintervalle (T = 25 undT = 30) kein Einpegeln der Differenzkak−ãkkzu sehen ist. Dies ist aber der Abhängigkeit der Fehlerkonstanten c=c(T) vom der Länge des Integrationsintervalls zuzuschreiben: Da sie mit wachsender Integrations- dauer wächst, vergrößert das auch den Gesamtfehler kak(T)−a(T)k=c(T) h^p+O(h^p+1) mita(T) als exakter Lösung zum ZeitpunktT. Weiter unten betrachten wir noch die Werte der Lösungskom- ponenten in Abhängigkeit vonT, dort sieht man das Konvergenz für grössereT kleinere Schrittweiten erfordert und mit dem Einpegeln der Lösung sich auch der Anstand kak−a˜kk zwischen ungestörter und gestörter Lösung stabilisiert.

Die nächste Tabelle zeigt die Auswirkung der Störung ∆ = 10⁻⁴ auf die Lösungak als Vielfaches der Störung ∆ = 10⁻⁴:

T = 20 T = 25 T = 30

k kak−˜akk/∆ kak−a˜kk/∆ kak−˜akk/∆

5 0.28 0.81 219.58

6 1.28 39.51 351.84

7 1.77 286.49 2629.05

8 0.97 2893.85 1955.72

9 1.21 1162.79 1538.78

10 1.44 78.18 1284.66

11 1.54 131.26 307.20

12 1.57 331.88 1330.85

13 1.58 967.22 1232.87

14 1.58 1320.21 1629.98

15 1.58 1426.76 1844.07

16 1.58 1454.58 1441.36

17 1.58 1461.61 1284.30

18 1.58 1463.37 1240.63

19 1.58 1463.82 1229.41

20 1.58 1463.93 1226.59

21 1.58 1463.95 1225.88

22 1.58 1463.96 1225.70

23 1.58 1463.96 1225.67

Man sieht leicht, das die Störung für große Integrationsintervalle eine sehr große Reaktion in der Lösung

˜

ak nach sich zieht. Hier zeigt sich das chaotische Verhalten des Lorenzsystems: Man kann von einer Trajektorie f¨ur einen Anfangswerta⁰nicht auf die Trajektorie zu einem leicht gest¨orten Anfangswert

˜

a⁰ schließen (zumindest für interessante, größere Zeiträume).

∆ = 10⁻⁶

(11)

T = 20 T = 25 T = 30 k ka_k−˜a_kk/∆ ka_k−a˜_kk/∆ ka_k−˜a_kk/∆

5 0.28 0.81 172.08

6 1.28 33.46 4537.51

7 1.76 357.23 9476.35

8 0.97 2135.91 5990.69

9 1.21 1410.67 14038.52

10 1.44 93.55 12841.71

11 1.54 112.16 8942.98

12 1.57 411.42 23801.39

13 1.58 1286.62 19700.81

14 1.58 1757.55 13364.06

15 1.59 1891.64 10601.62

∆ = 10⁻⁸

T = 20 T = 25 T = 30

k kak−˜akk/∆ kak−a˜kk/∆ kak−˜akk/∆

5 0.28 0.81 171.69

6 1.28 33.41 4485.52

7 1.76 358.13 10241.60

8 0.97 2127.72 7029.37

9 1.21 1412.73 13478.77

10 1.44 93.69 12266.99

11 1.54 111.99 9728.83

12 1.57 412.27 31742.48

13 1.58 1289.67 19501.67

14 1.59 1763.87 14269.47

15 1.59 1899.33 9525.28

An den Tabellen für ∆ = 10⁻⁶ und ∆ = 10⁻⁸ kann man erkennen, dass die Reaktion nicht von der Größe der Störung und der Schrittweitehk abhängt: Auch auf kleine Störungen in der Größenordnung von 10⁻⁸ reagiert die Lösung beiT = 30 mit einem Faktor von ca 10⁴.

Grafische Darstellung der Auswirkung einer St¨orung von 10⁻⁴in der ersten Komponente des Anfangs- wertes (0.1,0.1,0.1) auf die L¨osung:

-20 -15 -10 -5 0 5 10 15 20-25-20-15-10-5 0 5 10 15 20 25 5

10 15 20 25 30 35 40 45

nicht gestoert gestoert Startpunkt bei t=21.5 Loesung nicht gestoert bei t=28 Loesung gestoert bei t=28

(12)

In der Grafik ist die Entwicklung der Lösungen für den ursprünglichen und den mit ∆ = 10⁻⁴gestörten Startwert zu sehen. Die Lösungstrakjektorien unterscheiden sich fürt ∈[0,21.5] nur geringfügig (un- gefähr in der Größenordnung der Störung ∆ = 10⁻⁴). Deshalb wurden sie auch nicht in obiger Grafik dargestellt.

Das Kreuz in der Grafik stellt die beide Lösungen zum Zeitpunktt= 21.5 dar, die Trajektorien stimmen auch noch für einige Zeit überein (Schlaufe nach rechts oben). Mit fortschreitender Zeit separieren sich beide allerdings deutlich, zum Zeitpunkt t = 28 (markiert durch die Quadrate) befinden sich die Lösungen an völlig verschiedenen Positionen. (Im weiteren Verlauf nähern und entfernen sich auch wieder.)

Erstelle eine graphische Darstellung der L¨osung (3-dimensional!!) vom Startpunkt(0.1,0.1,0.1) f¨ur eink >7.

L¨osungstrajektorie f¨urt= 0..30 von Anfangswert (0.1,0.1,0.1):

-20 -15 -10 -5 0 5 10 15 20-30

-20 -10

0

10 20 30 0 5

10 15 20 25 30 35 40 45 50

’a3-30-9.res’

Nun untersuchen wir noch kurz den Einfluss der Intergrationsdauer auf die Konvergenzeigenschaften.

Dazu betrachten wir die Werte von x_k f¨urk= 5, . . . ,23 mit ∆ = 10⁻⁴:

T = 20 T = 25 T = 30

k xk x˜k xk x˜k xk x˜k

5 -3.494566 -3.494338 1.952341 1.947815 6.053674 7.104504 6 -11.270226 -11.266371 10.996204 10.798048 -6.480503 -5.254567 7 1.880002 1.878429 -12.676321 -12.640897 4.371818 -9.015788 8 -1.866830 -1.865680 -16.961988 -7.470797 12.283833 10.588341 9 -2.039448 -2.039303 8.830243 13.074962 -1.352231 6.396224 10 -2.071395 -2.072062 12.844207 12.428216 -1.121606 5.339958 11 -2.065532 -2.066553 12.544303 12.712908 -4.388492 -2.978470 12 -2.062175 -2.063306 15.836561 14.950599 -2.353931 0.515525 13 -2.061164 -2.062325 14.834763 16.322094 -4.621451 -9.523835 14 -2.060897 -2.062066 12.315842 15.856262 3.085576 -5.045878 15 -2.060829 -2.062000 11.408694 15.607651 -8.224703 1.156648 16 -2.060812 -2.061983 11.161890 15.534822 -8.444340 -1.205992 17 -2.060807 -2.061979 11.098822 15.515879 -8.158642 -1.722782 18 -2.060806 -2.061978 11.082958 15.511093 -8.057922 -1.844591 19 -2.060806 -2.061978 11.078985 15.509893 -8.030658 -1.874562 20 -2.060806 -2.061978 11.077991 15.509593 -8.023705 -1.882003 21 -2.060806 -2.061978 11.077742 15.509517 -8.021958 -1.883893 22 -2.060806 -2.061978 11.077680 15.509499 -8.021521 -1.884355 23 -2.060806 -2.061978 11.077665 15.509496 -8.021412 -1.884416

(13)

T = 20 T = 25 T = 30 k x_k x˜_k x_k x˜_k x_k x˜_k

Weiter oben haben wir f¨urT = 1 undk >5 beobachtet, dass die berechnete L¨osungakund insbesondere deren erste Komponentexk gegen einen Grenzwert konvergieren.

FürT = 20 sehen wir hier nun, dasxk sich fürk >9 stabilisiert und mit fortlaufendemk Konvergenz zeigt (die Komponentenyk undzk zeigen ähnliches Verhalten), d.h. die Zahl der sich für wachsendes k nicht mehr verändernden Dezimalstellen steigt (für gestörte und ungestörte Anfangswerte gleicher- maßen).

Im Falle vonT = 25 hat sich beik= 15 nur eine Dezimalstelle eingepegelt, die Rechnung mit kleineren Schrittweiten zeigt fortschreitenden Stabilisierung der numerischen L¨osung.

BeiT = 30 ist bei einer Schrittweitehk= 1/2¹⁵= 1/32768 =.000030517578125 noch nicht einmal das Vorzeichen der L¨osung erkennbar!! Beik= 23 sind nur die ersten 4 Dezimalstellen fixiert.

M¨ogliches C-Programm:

1 #include <math . h>

2

3 #define MINK 1 4 #define MAXK 2 0 5

6

7 typedef s t r u c t { double x , y , z ; } v a l ; 8

9 // r e c h t e S e i t e d e r ODE

10 double dxdt (double x , double y , double z ) { 11 return −10∗x + 1 0∗y ;

12 }

13 double dydt (double x , double y , double z ) { 14 return 28∗x − y − x∗z ;

15 }

16 double d z d t (double x , double y , double z ) { 17 return −8 . 0∗z / 3 . 0 + x∗y ;

18 } 19 20

21 double norm2 ( v a l a1 , v a l a0 ) { 22 double s ;

23 return s q r t ( ( a1 . x−a0 . x ) ∗ ( a1 . x−a0 . x ) 24 + ( a1 . y−a0 . y ) ∗ ( a1 . y−a0 . y ) 25 + ( a1 . z−a0 . z ) ∗ ( a1 . z−a0 . z ) ) ; 26 }

27

28 double e r r c ( v a l a1 , v a l a0 , i n t k ) {

29 return 4 ∗ norm2 ( a1 , a0 ) ∗ pow ( 2 , 2∗k ) / 3 ; 30 }

31 32 33

34 void k l o o p m i d p o i n t ( i n t kmin , i n t kmax ,

35 double T , double T0 ,

36 double X0 , double Y0 , double Z0 , 37 v a l∗ VT, i n t p r i n t )

38 {

39 long long i n t k , n , i , N;

40 double h , tk , xk , yk , zk , dxk , dyk , dzk ; 41

42 i f ( p r i n t )

43 p r i n t f ( ”## T % l f X0 % l f Y0 % l f Z0 % l f kmx %d\n” ,

44 T , X0 , Y0 , Z0 , kmax ) ;

45

(14)

46 n=1;

47 n = 1 < < (kmin−1 ) ; 48

49 f o r ( k = kmin ; k <= kmax ; k ++ ) {

50 n ∗= 2 ;

51 N = (long) c e i l ( f a b s (T−T0 )∗(double) n ) ; 52 h = (T−T0 ) / (double) N;

53

54 i f ( p r i n t ) p r i n t f ( ”#\n” ) ;

55 p r i n t f ( ”# k=%d n=%d h=%l f\n” , k , n , h ) ; 56 i f ( p r i n t ) p r i n t f ( ”#\n” ) ;

57

58 t k = T0 ;

59 xk = X0 ; yk = Y0 ; zk = Z0 ;

60

61 i f ( p r i n t ) p r i n t f ( ” % l f % l f % l f % l f\n” , xk , yk , zk , t k ) ; 62

63 f o r ( i = 1 ; i <N+ 1 ; i ++ ) { 64

65 // m i d p o i n t

66 dxk = xk + h∗dxdt ( xk , yk , zk ) / 2 . 0 ; 67 dyk = yk + h∗dydt ( xk , yk , zk ) / 2 . 0 ; 68 dzk = zk + h∗d z d t ( xk , yk , zk ) / 2 . 0 ; 69

70 xk += h ∗ dxdt ( dxk , dyk , dzk ) ;

71 yk += h ∗ dydt ( dxk , dyk , dzk ) ;

72 zk += h ∗ d z d t ( dxk , dyk , dzk ) ; 73

74 // t i m e s t e p

75 t k = i∗h ;

76

77 i f ( p r i n t ) p r i n t f ( ”% l f % l f % l f % l f\n” , xk , yk , zk , t k ) ;

78 }

79

80 VT[ k ] . x = xk ; VT[ k ] . y = yk ; VT[ k ] . z = zk ; 81

82 i f ( p r i n t ) p r i n t f ( ”\n\n” ) ;

83 e l s e {

84 p r i n t f ( ”% l f % l f % l f % l f ” , xk , yk , zk , t k ) ;

85 i f ( k > kmin ) p r i n t f ( ” %7.1 e\n” , norm2 (VT[ k ] ,VT[ k−1 ] ) ) ; 86 e l s e p r i n t f ( ”−−−−\n” ) ;

87 }

88 }

89 } 90

91 i n t main ( void ) { 92 double T = 1 . 0 ; 93 double T0 = 0 . 0 ; 94

95 double X0 = 0 . 1 ; 96 double Y0 = 0 . 1 ; 97 double Z0 = 0 . 1 ; 98

99 i n t i , p r i n t = 0 ; 100

101 v a l YT[MAXK+ 1 ] ; 102

103 p r i n t f ( ”## T % l f X0 % l f Y0 % l f ”

104 ”Z0 % l f kmx %d\n” ,

105 T , X0 , Y0 , Z0 , MAXK ) ;

106

(15)

107 k l o o p m i d p o i n t ( MINK, MAXK, T , T0 , X0 , Y0 , Z0 , YT, p r i n t ) ; 108

109 p r i n t f ( ”\n\n” ) ;

110 f o r ( i=MINK ; i <MAXK; i ++ ) 111 p r i n t f ( ” %4d % l f\n” , i , 112 e r r c (YT[ i ] ,YT[ i + 1 ] , i ) ) ; 113 }

phone: 030/2093-5820 fax: 030/2093-5859 e-mail: griewank@math.hu-berlin.de riehme@math.hu-berlin.de http://www.mathematik.hu-berlin.de/∼gaggle/MATHINF