WS13/14 NumerischeMethodenvongewöhnlichenDiﬀerentialgleichungen(AWP)Prof.Dr.-Ing.K.Warendorf,Prof.Dr.-Ing.P.Wolfsteiner NumerischeVerfahren

(1)

Explizite Zeitschrittverfahren für DGL 1. Ordnung Implizite Zeitschrittverfahren für DGL 1. Ordnung Anwendung auf DGL-Systeme 1. Ordnung Zustandsform

Numerische Verfahren

Numerische Methoden von gewöhnlichen Differentialgleichungen (AWP)

Prof. Dr.-Ing. K. Warendorf, Prof. Dr.-Ing. P. Wolfsteiner

Hochschule für Angewandte Wissenschaften München Fakultät 03

WS 13/14

(2)

Inhalt

1 Explizite Zeitschrittverfahren für DGL 1. Ordnung Einleitung

Explizites Euler-Verfahren

Runge-Kutta-Verfahren 4. Ordnung

2 Implizite Zeitschrittverfahren für DGL 1. Ordnung

3 Anwendung auf DGL-Systeme 1. Ordnung DGL-Systeme 1. Ordnung

4 Zustandsform

Reduktion von DGLn höherer Ordnung

Reduktion von Systemen von DGLn 2. Ordnung

(3)

Einleitung

Explizites Euler-Verfahren Runge-Kutta-Verfahren 4. Ordnung

Einleitung

Gegeben sei das Anfangswertproblem:

˙

y = f (t, y), y(t

0

) = y

0

Aufgabe

Berechnung eines Näherungswertes an der Stelle t

n

= b.

y

n

≈ y(b)

mit der Diskretisierung:

Zerlegung des Integrationsintervalles [a = t

0

, b] in: a = t

0

< t

1

< · · · < t

n

= b Feste Schrittweite: h =

^b−a_n

= t

i

− t

i−1

Achtung: Der Stabilitätsbereich eines expliziten Einschrittverfahrens und damit die Schrittweite h sind beschränkt. Für zu große h divergiert das Verfahren.

h ist auch nach unten beschränkt, da sonst die Rundungsfehler überwiegen!

(4)

Einleitung

Explizites Euler-Verfahren

Fortschreiterichtung: Steigung in t _i−1 Euler-Verfahren

Gegeben: Anfangswert: (t ₀ , y ₀ ), Schrittweite: h y _i = y _i−1 + h f (t _i−1 , y _i−1 ), i = 1, ..., n Globaler Diskretisierungsfehler: O(h).

⇒ Verfahren 1. Ordnung (Konvergenz)

(5)

Einleitung

Runge-Kutta-Verfahren 4. Ordnung

Fortschreiterichtung: mittlere Steigung: Berechnet aus den 4 Steigungen in t

i−1

, t

i

und 2 Punkten in der Mitte

Runge-Kutta-Verfahren 4. Ordnung

Gegeben: Anfangswert: (t

0

, y

0

), Schrittweite: h y

i

= y

i−1

+

¹₆

(k

1

+ 2k

2

+ 2k

3

+ k

4

), i = 1, ..., n k

1

= h f (t

i−1

, y

i−1

)

k

2

= h f (t

i−1

+

¹₂

h, y

i−1

+

¹₂

k

1

) k

3

= h f (t

i−1

+

¹₂

h, y

i−1

+

¹₂

k

2

) k

4

= h f (t

i−1

+ h, y

i−1

+ k

3

) Globaler Diskretisierungsfehler: O(h

⁴

).

⇒ Verfahren 4. Ordnung (Konvergenz)

Matlab: ode45(...). Dieses und entsprechende Verfahren sind mit automatischer

Schrittweitensteuerung implementiert.

(6)

Implizites Euler-Verfahren

Implizites Verfahren: Die Funktion f wird nicht an dem schon berechneten Punkt (t

i−1

, y

i−1

) ausgewertet, sondern an dem neuen Punkt (t

i

, y

i

)

Implizites Euler-Verfahren

Gegeben: Anfangswert: (t

0

, y

0

), Schrittweite: h y

i

= y

i−1

+ h f (t

i

, y

i

), i = 1, ..., n Globaler Diskretisierungsfehler: O(h).

⇒ Verfahren 1. Ordnung (Konvergenz) Das Verfahren ist unbeschränkt absolut stabil (A-stabil), d.h. h unterliegt aufgrund der Stabilität keiner Beschränkung nach oben.

Es entsteht eine (nicht)-lineare Gleichung (bzw. System). Diese muss mit geeigneten (numerischen) Verfahren (z.B.) Newton-Verfahren, Fixpunktiteration) gelöst werden.

Implizite Verfahren werden insbesondere zur Lösung von steifen DGL-Systemen genutzt.

(7)

DGL-Systeme 1. Ordnung

DGL-Systeme 1. Ordnung

Die numerischen Verfahren zur Lösung von DGLn 1. Ordnung lassen sich entsprechend auf Systeme 1. Ordnung übertragen. Die skalaren

Operationen müssen als Vektoroperationen aufgefasst werden.

DGL-System 1. Ordnung

˙

y ₁ = f ₁ (t, y ₁ , . . . , y _n )

˙

y ₂ = f ₂ (t, y ₁ , . . . , y _n ) .. .

˙

y _n = f _n (t , y ₁ , . . . , y _n )



 



 



: Vektorschreibweise: ~ y ˙ = F(t , ~ y)

mit den Anfangsbedingungen

y ₁ (t ₀ ) = y ₁ ⁽⁰⁾ , y ₂ (t ₀ ) = y ₂ ⁽⁰⁾ , . . . , y _n (t ₀ ) = y _n ⁽⁰⁾ bzw. ~ y(t ₀ ) = ~ y ₀ .

(8)

DGL-Systeme 1. Ordnung

Explizites Euler-Verfahren für DGL-System 1. Ordnung

Euler-Verfahren

Gegeben: Anfangswerte: y ₁ (t ₀ ) = y ⁽⁰⁾ ₁ , y ₂ (t ₀ ) = y ₂ ⁽⁰⁾ , . . . , y _n (t ₀ ) = y _n ⁽⁰⁾ bzw. ~ y(t ₀ ) = ~ y ₀ , Schrittweite: h

~ y _i = ~ y _i−1 + h · F(t _i−1 , ~ y _i−1 ), i = 1, ..., n

(9)

Reduktion von DGLn höherer Ordnung Reduktion von Systemen von DGLn 2. Ordnung Reduktion von Systemen von DGLn 2. Ordnung

Reduktion von DGLn höherer Ordnung auf ein System 1. Ordnung (Zustandsform) Die numerische Lösung von DGLn höherer Ordnung oder von Systemen, erfordert ein Umschreiben in die sogenannte Zustandsform (System 1. Ordnung).

DGL n-ter Ordnung:

y

⁽ⁿ⁾

= f (t, y, y, . . . , ˙ y

⁽ⁿ⁻¹⁾

), Anfangsbedingungen: y(t

0

), y(t ˙

0

), . . . , y

⁽ⁿ⁻¹⁾

(t

0

) Zustandsform

Zustandsgrößen: z

1

= y, z

2

= ˙ y, z

3

= ¨ y, . . . z

n

= y

⁽ⁿ⁻¹⁾

Umformung in ein System 1. Ordnung:

˙

z

1

= z

2

, z

1

(t

0

) = y(t

0

)

˙

z

2

= z

3

, z

2

(t

0

) = ˙ y(t

0

) . .

. , . . .

˙

z

n−1

= z

n

, z

n−1

(t

0

) = y

⁽ⁿ⁻²⁾

(t

0

)

˙

z

n

= f (t, z

1

, . . . , z

n

) , z

n

(t

0

) = y

⁽ⁿ⁻¹⁾

(t

0

)

(10)

Zustandsform für DGLn 2. Ordnung I

DGL 2. Ordnung:

¨

y = f (t, y, y) ˙ mit AB: y(t ₀ ) = y ₁ ⁽⁰⁾ , y(t ˙ ₀ ) = y ⁽⁰⁾ ₂

lässt sich durch z ₁ = y, z ₂ = ˙ y überführen in die Zustandsform:

˙

z ₁ = z ₂ , z ₁ (t ₀ ) = z ₁ ⁽⁰⁾ = y ₁ ⁽⁰⁾

˙

z ₂ = f (t, z ₁ , z ₂ ) , z ₂ (t ₀ ) = z ₂ ⁽⁰⁾ = y ₂ ⁽⁰⁾ (1)

Vektorielle Schreibweise: ~ z ˙ =



 z ₂ f (t,~ z)





(11)

Zustandsform für DGLn 2. Ordnung II

Euler-Verfahren für obiges System (1)

Gegeben: Anfangswerte: (t ₀ , z ₁ ⁽⁰⁾ ),(t ₀ , z ₂ ⁽⁰⁾ ) Schrittweite: h , Schrittanzahl: m

z ₁ ⁽ⁱ⁾ = z ₁ ⁽ⁱ⁻¹⁾ + h · z ₂ ⁽ⁱ⁻¹⁾

z ₂ ⁽ⁱ⁾ = z ₂ ⁽ⁱ⁻¹⁾ + h · f (t i−1 , z ₁ ⁽ⁱ⁻¹⁾ , z ₂ ⁽ⁱ⁻¹⁾ ) i = 1, ..., m

Vektorielle Schreibweise: ~ z _i = ~ z _i−1 + h · ~ z ˙ _i−1

(12)

Zustandsform einer DGL 2. Ordnung: Punktpendel

Nichtlineare DGL 2. Ordnung in ϕ:

m · l

²

· ϕ(t) = ¨ −m · g · l · sin (ϕ(t))

⇒ ϕ(t) = ¨ − g

l · sin (ϕ(t)) Anfangsbedingungen:

Anfangswinkel: ϕ(0), Anfangswinkelgeschwindigkeit: ϕ(0) ˙ Zustandsform

Zustandsgrößen:

Lage: z

1

(t) = ϕ(t), Geschwindigkeit: z

2

(t) = ˙ ϕ(t) Umformung in ein System 1. Ordnung:

˙

z

1

= z

2

, z

1

(0) = ϕ(0)

˙

z

2

= −

^g_l

· sin (z

1

(t)) , z

2

(0) = ˙ ϕ(0)

(13)

Zustandsform eines Systems 2. Ordnung: I

System von DGL 2. Ordnung in y

1

, . . . , y

n

:

¨

y

1

= f

1

(t, y

1

, . . . , y

n

, y ˙

1

. . . , y ˙

n

) . .

.

¨

y

n

= f

n

(t, y

1

, . . . , y

n

, y ˙

1

. . . , y ˙

n

) Anfangsbedingungen:

y

1

(t

0

), . . . , y

n

(t

0

),

˙

y

1

(t

0

), . . . , y ˙

n

(t

0

)

(14)

Zustandsform eines Systems 2. Ordnung: II

Zustandsform Zustandsgrößen:

z

1

= y

1

, . . . , z

n

= y

n

z

n+1

= ˙ y

1

, . . . , z

2n

= ˙ y

n

Umformung in ein System 1. Ordnung:

˙

z

1

= z

n+1

, z

1

(t

0

) = y(t

0

)

. .

. , . . .

˙

z

n

= z

2n

, z

n

(t

0

) = y

n

(t

0

)

˙

z

n+1

= f

1

(t, z

1

, . . . , z

n

, z

n+1

. . . , z

2n

) , z

n+1

(t

0

) = ˙ y

1

(t

0

) . .

. , . . .

(15)

Zustandsform eines Systems 2. Ordnung: Federpendel I

System von DGL 2. Ordnung in ϕ, l:

m · (l(t) ¨ ϕ(t) + 2˙ l(t) ˙ ϕ(t)) = −mg sin(ϕ(t))

m · (¨ l(t) − l(t) ˙ ϕ

²

(t)) = −c(l(t) − l

F

) + mg cos(ϕ(t))

⇒

¨

ϕ(t) = 1

l (−g sin(ϕ(t)) − 2˙ l(t) ˙ ϕ(t))

¨ l(t) = − c

m (l(t) − l

F

) + g cos(ϕ(t)) + l(t) ˙ ϕ

²

(t) Anfangsbedingungen:

Anfangslage: ϕ(0), l(0)

Anfangswinkelgeschwindigkeit: ϕ(0), ˙ ˙ l(0)

(16)

Zustandsform eines Systems 2. Ordnung: Federpendel II

Zustandsform Zustandsgrößen:

Lage: z

1

(t) = ϕ(t), z

2

(t) = l(t) Geschwindigkeit: z

3

(t) = ˙ ϕ(t), z

4

(t) = ˙ l(t) Umformung in ein System 1. Ordnung:

˙

z

1

= z

3

, z

1

(0) = ϕ(0)

˙

z

2

= z

4

, z

2

(0) = l(0)

˙

z

3

= −

_z¹

2

(g sin(z

1

) + 2z

4

z

3

) , z

3

(0) = ˙ ϕ(0)

˙

z

4

= −

_m^c

(z

2

− l

F

) + g cos(z

1

) + z

2

z

₃²

, z

4

(0) = ˙ l(0)

(17)

Zustandsform eines Systems 2. Ordnung: Massenmatrix I

Falls das DGL-System nicht in expliziter Form gegeben ist und sich nur schwer nach den 2. Ableitungen auflösen lässt (starke Kopplung), ist eine andere Vorgehensweise zu empfehlen.

System von 2 DGL (linear in den 2. Ableitungen) a

11

¨ y

1

+ a

12

¨ y

2

= f

1

(t, y

1

, y

2

, y ˙

1

, y ˙

2

) a

21

¨ y

1

+ a

22

¨ y

2

= f

2

(t, y

1

, y

2

, y ˙

1

, y ˙

2

)

a

ij

können von t, y

1

, y

2

, y ˙

1

und y ˙

2

aber nicht von ¨ y

1

und ¨ y

2

abhängen.

(18)

Zustandsform eines Systems 2. Ordnung: Massenmatrix II

Zustandsform

Zustandsgrößen: z

1

= y

1

,z

2

= y

2

, z

3

= ˙ y

1

,z

4

= ˙ y

2

˙

z

1

= z

3

, z

1

(t

0

) = y

1

(t

0

)

˙

z

2

= z

4

, z

2

(t

0

) = y

2

(t

0

)

a

11

z ˙

3

+ a

12

z ˙

4

= f

1

(t,z

1

,z

2

,z

3

, z

4

) , z

3

(t

0

) = ˙ y

1

(t

0

) a

21

z ˙

3

+ a

22

z ˙

4

= f

2

(t,z

1

,z

2

,z

3

, z

4

) , z

4

(t

0

) = ˙ y

2

(t

0

) Matrix-Schreibweise:







1 0 0 0

0 1 0 0

0 0 a

11

a

12

0 0 a

21

a

22







| {z }

Massenmatrix:M=M(t,~z)

· 





˙ z

1

˙ z

2

˙ z

3

˙ z

4







| {z }

~z˙

=







z

3

z

4

f

1

(t,z

1

,z

2

,z

3

, z

4

) f

2

(t,z

1

,z

2

,z

3

, z

4

)







| {z }

F(t, ~z)

Matlab kann diese gekoppelten Differentialgleichungen lösen. Dazu muss M mit dem Befehl

odeset(’Mass’, @Matrix) bekannt gemacht werden.

WS13/14 NumerischeMethodenvongewöhnlichenDiﬀerentialgleichungen(AWP)Prof.Dr.-Ing.K.Warendorf,Prof.Dr.-Ing.P.Wolfsteiner NumerischeVerfahren

Numerische Verfahren

Numerische Methoden von gewöhnlichen Differentialgleichungen (AWP)

Prof. Dr.-Ing. K. Warendorf, Prof. Dr.-Ing. P. Wolfsteiner

Hochschule für Angewandte Wissenschaften München Fakultät 03

WS 13/14

Inhalt

1 Explizite Zeitschrittverfahren für DGL 1. Ordnung Einleitung

Explizites Euler-Verfahren

Runge-Kutta-Verfahren 4. Ordnung

2 Implizite Zeitschrittverfahren für DGL 1. Ordnung

3 Anwendung auf DGL-Systeme 1. Ordnung DGL-Systeme 1. Ordnung

4 Zustandsform

Reduktion von DGLn höherer Ordnung

Reduktion von Systemen von DGLn 2. Ordnung

Reduktion von Systemen von DGLn 2. Ordnung

Einleitung

Gegeben sei das Anfangswertproblem:

˙

y = f (t, y), y(t

) = y

Aufgabe

Berechnung eines Näherungswertes an der Stelle t

= b.

y

≈ y(b)

mit der Diskretisierung:

Zerlegung des Integrationsintervalles [a = t

, b] in: a = t

< t

< · · · < t

= b Feste Schrittweite: h =

= t

− t

Achtung: Der Stabilitätsbereich eines expliziten Einschrittverfahrens und damit die Schrittweite h sind beschränkt. Für zu große h divergiert das Verfahren.

h ist auch nach unten beschränkt, da sonst die Rundungsfehler überwiegen!

Explizites Euler-Verfahren

Fortschreiterichtung: Steigung in t i−1 Euler-Verfahren

Gegeben: Anfangswert: (t 0 , y 0 ), Schrittweite: h y i = y i−1 + h f (t i−1 , y i−1 ), i = 1, ..., n Globaler Diskretisierungsfehler: O(h).

⇒ Verfahren 1. Ordnung (Konvergenz)

Runge-Kutta-Verfahren 4. Ordnung

Fortschreiterichtung: mittlere Steigung: Berechnet aus den 4 Steigungen in t

, t

und 2 Punkten in der Mitte

Runge-Kutta-Verfahren 4. Ordnung

Gegeben: Anfangswert: (t

, y

), Schrittweite: h y

= y

+

(k

+ 2k

+ 2k

+ k

), i = 1, ..., n k

= h f (t

, y

)

k

= h f (t

+

h, y

+

k

) k

= h f (t

+

h, y

+

k

) k

= h f (t

+ h, y

+ k

) Globaler Diskretisierungsfehler: O(h

).

⇒ Verfahren 4. Ordnung (Konvergenz)

Matlab: ode45(...). Dieses und entsprechende Verfahren sind mit automatischer

Schrittweitensteuerung implementiert.

Implizites Euler-Verfahren

Fortschreiterichtung: Steigung in t _i−1 Euler-Verfahren

Gegeben: Anfangswert: (t ₀ , y ₀ ), Schrittweite: h y _i = y _i−1 + h f (t _i−1 , y _i−1 ), i = 1, ..., n Globaler Diskretisierungsfehler: O(h).

y ₁ = f ₁ (t, y ₁ , . . . , y _n )

y ₂ = f ₂ (t, y ₁ , . . . , y _n ) .. .

y _n = f _n (t , y ₁ , . . . , y _n )

y ₁ (t ₀ ) = y ₁ ⁽⁰⁾ , y ₂ (t ₀ ) = y ₂ ⁽⁰⁾ , . . . , y _n (t ₀ ) = y _n ⁽⁰⁾ bzw. ~ y(t ₀ ) = ~ y ₀ .

Gegeben: Anfangswerte: y ₁ (t ₀ ) = y ⁽⁰⁾ ₁ , y ₂ (t ₀ ) = y ₂ ⁽⁰⁾ , . . . , y _n (t ₀ ) = y _n ⁽⁰⁾ bzw. ~ y(t ₀ ) = ~ y ₀ , Schrittweite: h

~ y _i = ~ y _i−1 + h · F(t _i−1 , ~ y _i−1 ), i = 1, ..., n