12 13

Numerik f¨ ur Ingenieure und Naturwissenschaftler

Wolfgang Dahmen und Arnold Reusken

ISBN 3-540-25544-3, Springer Verlag

Inhaltsangabe 1. Einleitung

2. Fehleranalyse: Kondition, Rundungsfehler, Stabilit¨at 3. Lineare Gleichungssysteme, direkte L¨osungsverfahren 4. Lineare Ausgleichsrechnung

5. Nichtlineare Gleichungssysteme, iterative L¨osungsverfahren 6. Nichtlineare Ausgleichsrechnung

7. Berechnung von Eigenwerten 8. Interpolation

9. Splinefunktionen

10. Numerische Integration

11. Gew¨ohnliche Differentialgleichungen 12. Partielle Differentialgleichungen

13. Große d¨unnbesetzte lineare Gleichungssysteme

14. Numerische Simulationen: Vom Pendel bis zum Airbus

Ubersicht¨

2

12 13

11 10 8 5 7

3 6

4

9

5.7 4.5-4.7

11.6.2, 11.8.5, 11.9.3

12.2, 12.4

8.2.7, 8.4, 8.5

Korrekturen

• Seite 12, Zeile -2: y statt u.

• Seite 31, Bemerkung 2.26 (ii): k · kX, k · kY statt k · kV, k · kW.

• Seite 32, Zeile -3: kbk∞ statt kxk∞.

• Seite 76: (3.56)

L˜⁻¹₁ · · ·L˜⁻¹_n−1R statt L˜⁻¹₁ · · ·L˜_n−1R

• Seite 161, Zeile 3: (3.13) statt 3.13.

• Seite 183, Zeile -12: (5.43) statt 5.43.

• Seite 271, Zeile -8: (8.11) statt 8.11.

• Seite 282, Zeilen -2 und -8: erg¨anze “f¨ur alle x ∈ [0,1]”.

• Seite 291, Zeile 12: (8.52) statt (8.27).

• Seite 298, Zeile -1: Schrittweite statt Schrittweise.

Beispiel 1.2.

”Technische Aufgabe“ :

Konstruktion eines Taktmechanismus zu einer vorgegebenen Taktzeit T.

Ansatz: Modellierung mit einem Pendel.

Aufgabe:

Bestimmung der erforderlichen Anfangsauslenkung zu einer vorgege- benen Schwingungsdauer T.

Die Dynamik des Pendels:

φ^′′(t) = −csin(φ(t)), c := g ℓ, mit Anfangsbedingungen

φ(0) = x, φ^′(0) = 0.

Parameter g, ℓ, x: Fallbeschleunigung, Pendell¨ange (ℓ = 0.6 m), Anfangsauslenkung x.

SS SS

SS SS

SS SS

SS

φ

?

mg

=

mg sinφ

SS SS S

SS SSS

Abbildung 1.3.: Bewegung f¨ur x = ^π₂.

−1.5

−1

−0.5 0 0.5 1 1.5 2

phi

Aufgabe: bestimme x = x^∗ ∈ (0, ^π₂), wof¨ur T(x^∗) = 1.8 gilt.

Das Pendel steht zum Zeitpunkt T /4 gerade senkrecht:

φ(T /4, x^∗) = 0.

Sei

f(x) := φ(T /4, x) = φ(0.45, x).

Problem der Takterkonstruktion: bestimme die Nullstelle x^∗ ∈ (0, ^π₂) dieser (nur implizit gegebenen) Funktion f:

f(x^∗) = 0.

Dies nennt man ein Nullstellenproblem.

Auswertung von f: L¨osung einer Anfangswertaufgabe

→ Diskretisierungsfehler.

Implementierung eines L¨osungsverfahrens → Rundungsfehler.

Weitere Fehler: Modellfehler, Datenfehler.

Abbildung 1.4.: Werte φ(0.45, x)˜ ≈ f(x)

−0.1

−0.08

−0.06

−0.04

−0.02 0 0.02 0.04 0.06

Ziele der Vorlesung

F¨ur unterschiedliche Problemstellungen (L¨osen eines linearen Glei- chungssystems, Berechnung eines Integrals, L¨osen einer Differential- gleichung) werden folgende Themen behandelt:

• Kondition (= Empfindlichkeit f¨ur St¨orungen) eines Problems.

• Wichtige numerische L¨osungsverfahren.

• Stabilit¨at (= Empfindlichkeit f¨ur St¨orungen) der L¨osungsverfah- ren.

• Effizienz (= Anzahl der Rechenoperationen, Speicherbedarf ) der L¨osungsverfahren, d.h., der numerische Aufwand, der n¨otig ist, um eine gew¨unschte

”L¨osungsqualit¨at“ , sprich Genauigkeit zu erzielen.

KAPITEL 10. Numerische Integration

10.1 Einleitung Sei

I =

Z _b

a f(x)dx, I˜=

Z _b a

f˜(x) dx.

Es gilt

|I − I˜| =

Z _b

a f(x) − f˜(x) dx ≤

Z _b

a |f(x) − f˜(x)|dx ≤ (b − a)kf − f˜k_∞.

|I − I˜|

≤ (b − a) kf − f˜k_∞

b =

R_b

a kfk_∞ dx

b · kf − f˜k_∞

=: κrel

kf − f˜k_∞ .

Die g¨angige Strategie zur n¨aherungsweisen Berechnung von

Z _b

a f(x) dx l¨aßt sich folgendermaßen umreißen:

1. Man unterteile [a, b] in Teilintervalle [t_k−1, t_k] z.B. mit t_j = a + jh, j = 0, . . . , n, h = ^b−_n^a.

2. Approximiere f auf jedem Intervall [t_k−1, t_k] durch eine einfach zu integrierende Funktion g_k, und verwende

n X k=1

Z _tk

t_k−1 g_k(x) dx ≈

n X k=1

Z _tk

t_k−1 f(x)dx =

Z _b

a f(x) dx als N¨aherung f¨ur das exakte Integral.

Trapezregel Dabei w¨ahlt man speziell

g_k(x) = x − t_k−1

h f(t_k) + t_k − x

h f(t_k−1),

d.h. die lineare Interpolation an den Intervallenden von [t_k−1, t_k].

Folglich ist ^R_t^tk

k−1 g_k(x) dx gerade die Fl¨ache h

2[f(t_k−1) + f(t_k)] .

s

s f

←− ^h₂[f(tk−1) + f(tk)]

······

············

··············

···············

················

·················

··················

···················

···············

··········

·····

Dies liefert die

summierte Trapezregel T(h) = h

1

2f(a) + f(t₁) + · · · + f(t_n−1) + 1

2f(b)

als N¨aherung f¨ur ^R_a^b f(x) dx.

F¨ur den Verfahrensfehler der Teilintegrale gilt folgende Darstellung:

Lemma 10.1. Sei f ∈ C²([t_k−1, t_k]). Es gilt:

h

2[f(t_k−1) + f(t_k)] =

Z _tk

t_k−1 f(x) dx + f^′′(ξ_k)

12 h³ , ξ_k ∈ [t_k−1, t_k].

F¨ur den Verfahrensfehler von T(h) ergibt sich damit die Absch¨atzung

T(h) −

Z _b

a f(x) dx =

n X k=1

f^′′(ξ_k) 12 h³

≤ h³ 12

n X k=1

f^′′(ξ_k) ≤ h³

12n max

x∈[a,b]

f^′′(x) . Mit nh = b − a ergibt sich insgesamt die Fehlerschranke

T(h) −

Z _b

a f(x) dx ≤ h²

12(b − a) max

x∈[a,b]

f^′′(x) .

Ebenfalls erh¨alt man wegen E(h) := T(h) −

Z _b

a f(x) dx =

n X k=1

f^′′(ξ_k)

12 h³ = h² 12

n X k=1

hf^′′(ξ_k) und

h→lim0

E(h)

h² = 1 12

Z _b

a f^′′(x)dx = 1 12

f^′(b) − f^′(a)

die Fehlersch¨atzung

E(h) ≈ Eˆ(h) := h² 12

f^′(b) − f^′(a) .

Beispiel 10.2.

Zur n¨aherungsweisen Berechnung von I =

Z _π/2

0 xcos x + e^x dx = π

2 + e^12π − 2

mit der Trapezregel ergeben sich die in folgender Tabelle angegebenen N¨aherungswerte, Verfahrensfehler und Fehlersch¨atzungen.

n T(h) |E(h)| = |T(h) − I| |E(h)ˆ | = ₁₂^h2|f^{′ π}₂ − f^′(0)|

4 4.396928 1.57e−02 1.59e−02

8 4.385239 3.97e−03 3.98e−03

16 4.382268 9.95e−04 9.96e−04

32 4.381523 2.49e−04 2.49e−04

10.2 Newton-Cotes-Formeln

F¨ur ein typisches Teilintervall [t_k−1, t_k] stehe der Einfachheit halber im Folgenden [c, d]. Seien nun

x₀, . . . , x_m ∈ [c, d] verschiedene Punkte.

Integration des Interpolationspolynoms liefert die Quadraturformel I_m(f) =

Z _d

c P(f|x₀, . . . , x_m)(x) dx.

Satz 10.3. Sei I_m(f) wie oben. F¨ur jedes Polynom Q ∈ Π_m gilt I_m(Q) =

Z _d

c Q(x) dx.

Man sagt, die Quadraturformel ist exakt vom Grade m.

Lemma 10.4. Es gibt Gewichte c₀, . . . , c_m, so daß I_m(f) die Form I_m(f) = h

m X j=0

c_jf(x_j) hat, wobei wieder h = d − c. Die c_j sind durch

c_j = 1 h

Z _d c

m Y

k=0 k6=j

x − x_k

x_j − x_k dx = 1 h

Z _d

c ℓ_jm(x)dx

gegeben, wobei ℓ_jm die Lagrange-Fundamentalpolynome zu den St¨utzstellen x₀, . . . , x_m sind.

W¨ahlt man speziell die St¨utzstellen x_j ¨aquidistant x₀ = c + 1

2h =: c + ξ₀h, wenn m = 0 ,

Man kann dann die Quadraturformel in der Form I_m(f) = h

m X j=0

c_jf(c + ξ_jh)

mit normierten St¨utzstellen ξ_j und Gewichten c_j schreiben, die jetzt unabh¨angig vom speziellen Intervall [c, d] sind. G¨angige Beispiele:

m ξ_j c_j I_m(f) − Rd

c f(x)dx 0 Mittelpunktsregel ¹

2 1 −₂₄¹ h³f⁽²⁾(ξ)

1 Trapezregel 0, 1 ¹

2, ¹

2

1

12h³f⁽²⁾(ξ) 2 Simpson-Regel 0, ¹₂,1 ¹

6, ⁴₆, ¹_{6 90}¹ (¹

2h)⁵f⁽⁴⁾(ξ)

3 ³

8-Regel 0, ¹₃, ²₃,1 ¹

8, ³₈, ³₈, ¹_{8 80}³ (¹

3h)⁵f⁽⁴⁾(ξ) 4 Milne-Regel 0, ¹

4, ¹

2, ³

4,1 ⁷

90, ³²

90, ¹²

90, ³²

90, ⁷

90

8 945(¹

4h)⁷f⁽⁶⁾(ξ)

Summierte Newton-Cotes-Formeln

Als Beispiel behandeln wir die summierte Simpson-Regel.

S(h) =

Z _b

a f(x) dx + E(h) mit

S(h) = h 6

f(t₀) + 4f

t₀ + t₁ 2

+ 2f(t₁) + 4f

t₁ + t₂ 2

+ 2f(t₂) + . . . + 2f(t_n−1) + 4f

t_n−1 + t_n 2

+ f(t_n)

und

n

X 1 1 _{5 (4)} h⁴ Xⁿ ₍₄₎

Es gilt, wegen nh = b − a,

|E(h)| ≤ h⁴

2880(b − a)kf⁽⁴⁾k_∞ , E(h) ≈ h⁴

2880

Z _b

a f⁽⁴⁾(x) dx = h⁴ 2880

f⁽³⁾(b) − f⁽³⁾(a).

Man beachte, daß beim Aufsummieren der einzelnen Teilintegrale,

Z _b

a f(x)dx =

n X k=1

Z _tk

t_k−1 f(x) dx, im Fehler eine h-Potenz verloren geht.

Beispiel 10.5

F¨ur das Integral in Beispiel 10.2 ergeben sich die Resultate wie in fol- gender Tabelle.

n S(h) |E(h)| ₂₈₈₀^h4 |f⁽³⁾(^π₂) − f⁽³⁾(0)| 4 4.381343022 6.93e−05 6.92e−05

8 4.381278035 4.33e−06 4.33e−06 16 4.381273978 2.70e−07 2.70e−07 32 4.381273725 1.69e−08 1.69e−08

10.3 Gauß-Quadratur

Zielvorgaben:

Entwickle f¨ur m ∈ N eine Formel

m X i=0

w_if(x_i) =

Z _d

c P(f|x₀, . . . , x_m)(x)dx mit:

• positiven Gewichten w_i, i = 0, . . . , m;

• mit m¨oglichst hohem Exaktheitsgrad n ≥ m, d.h.,

Z _d

c Q(x) dx =

m X i=0

w_iQ(x_i), ∀ Q ∈ Π_n.

Der Exaktheitsgrad bei Newton-Cotes-Formeln I_m(f) ist entweder m oder m + 1. Es zeigt sich, daß man dies verbessern kann.

Allerdings sieht man leicht, daß man mit einer Formel dieses Typs h¨ochstens den Exaktheitsgrad 2m + 1 realisieren kann.

Daß man jedoch den maximalen Exaktheitsgrad n = 2m+ 1 realisieren kann, zeigen die Gaußschen Quadraturformeln.

Satz 10.6. Sei m ≥ 0. Es existieren St¨utzstellen x₀, . . . , x_m ∈ (c, d) und positive Gewichte w₀, . . . , w_m, so daß mit h = d − c

h

m X i=0

w_if(x_i) =

Z _d

c f(x) dx + E_f(h) und

E_Q = 0 f¨ur alle Q ∈ Π_2m+1. Ferner gilt f¨ur passendes ξ ∈ [c, d]

4

Daß die Gewichte w_j tats¨achlich positiv sind, ergibt sich aus der Exaktheit vom Grade 2m + 1 durch Anwendung auf das spezielle Polynom

q(x) :=

m Y i=0,i6=k

(x − x_i)² ∈ Π_2m ⊂ Π_2m+1, denn

0 <

Z _d

c q(x)dx =

m X i=0

w_iq(x_i) = w_kq(x_k) = w_k

m Y i=0,i6=k

(x_k − x_i)².

Numerische Tests

Wir untersuchen nun den in der Fehlerformel auftretenden Faktor

C_k,h := (k!)⁴

((2k)!)³(2k + 1)h^2k+1

(k = m + 1). F¨ur glatte Funktionen (d.h., |f^(2k)| wird nicht allzu groß, wenn k gr¨oßer wird) wird die Qualit¨at der Gauß-Quadratur im Wesentlichen durch den Faktor C_k,h bestimmt.

h k = 2 k = 4 k = 8

4 2.4e−01 1.5e−04 2.9e−13 2 7.4e−03 2.9e−07 2.2e−18 1 2.3e−04 5.6e−10 1.7e−23

Sei

I_k,n ≈

Z _b

a f(x)dx = I(f)

die Quadraturformel, wobei [a, b] in n Teilintervalle mit L¨ange ^b−_n^a = h unterteilt wird und auf jedem Teilintervall eine Gauß-Quadratur mit k St¨utzstellen angewandt wird.

Sowohl f¨ur I_2k,n als auch f¨ur I_k,2n wird die Anzahl der Funktionsaus- wertungen etwa verdoppelt im Vergleich zu I_k,n.

In obiger Tabelle kann man sehen, daß man

I − I_2k,n ≪ I − I_k,2n erwarten darf.

Daher wird in der Praxis bei Gauß-Quadratur n in der Regel klein gew¨ahlt, oft sogar n = 1.

Beispiel 10.7.

Die Gauß-Quadratur mit [c, d] = [0, ^π₂] (d.h. n = 1) f¨ur das Integral in Beispiel 10.2 ergibt die Resultate:

m I_m |Im − I|

1 4.3690643196 1.22e−03 2 4.3813023502 2.86e−05 3 4.3812734352 2.73e−07 4 4.3812737083 5.18e−10

Man sieht, daß in diesem Beispiel die Genauigkeit der Gauß-Quadratur mit 5 Funktionswerten (m = 4; k = 5) besser ist als die der Simpson- Regel angewandt auf n = 32 Teilintervalle (vgl. Tabelle 10.3), wobei

Beispiel 10.8.

Es sei [c, d] = [−1,1] und m = 1. Die Gauß-Quadraturformel I₁(f) = 2(c₀f(x₀) + c₁f(x₁))

muß f¨ur p ∈ Π₃ exakt sein, d.h.

Z ₁

−1 p(x)dx = 2(c₀p(x₀) + c₁p(x₁)) f¨ur p(x) = x^k, k = 0,1,2,3.

Aus

Z ₁

−1 x^k dx = 2(c₀x^k₀ + c₁x^k₁), k = 0,1,2,3, erh¨alt man die Gleichungen

2 = 2(c₀ + c₁) , 0 = 2(c₀x₀ + c₁x₁) ,

2

3 = 2(c₀x²₀ + c₁x²₁) , 0 = 2(c₀x³₀ + c₁x³₁) .

Dieses nichtlineare Gleichungssystem hat genau zwei L¨osungen:

c₀ = c₁ = 1

2, x₀ = −1 3

√3, x₁ = 1 3

√3, c₀ = c₁ = 1

2, x₀ = 1 3

√3, x₁ = −1 3

√3.

10.4 Extrapolation und Romberg-Quadratur Zu berechnen sei das Integral

I =

Z _b

a f(x) dx.

Die Trapezsumme T(h) liefert eine Approximation der Ordnung h². Die wesentliche Grundlage f¨ur den Erfolg von Extrapolationstechni- ken bildet eine sogenannte asymptotische Entwicklung des Diskreti- sierungsfehlers:

T(h) − I = c₁h² + c₂h⁴ + c₃h⁶ + · · · + c_ph^2p + O(h^2p+2) . Hieraus erh¨alt man

T1

2h − I = c₁1

4h² + ˆc₂h⁴ + . . . + ˆc_ph^2p + O(h^2p+2) , also

Man kann diese Idee systematisch weitertreiben. Sei T₁(h) = 4T(¹₂h) − T(h)

3 .

Es gilt

T₁(h) − I = ˜c₁h⁴ + ˜c₂h⁶ + . . . + O(h^2p+2), und damit

T₁(1

2h) − I = ˜c₁ 1

16h⁴ + ˜c₂ 1

64h⁶ + . . . + O(h^2p+2), und

16 15

T₁(1

2h) − I − 1 15

T₁(h) − I = 16T₁(¹₂h) − T₁(h)

15 − I

= d₁h⁶ + d₂h⁸ + . . . + O(h^2p+2) . Man erkennt, daß die Entwicklung des Fehlers der Quadraturformel

T₂(h) := 16T₁(¹₂h) − T₁(h) 15

mit einem Glied der Ordnung h⁶ beginnt.

Allgemeine Idee der Extrapolation

Die zugrunde liegende Idee ist folgende. Der gesuchte Wert ist T(0) =

Z _b

a f(x)dx = I.

Bestimmt man das lineare Interpolationspolynom der Funktion x → g(x) = T(√

x) = I + c₁x + c₂x² + . . . + c_px^p + O(x^p+1), (x ↓ 0),

zu den Punkten (h², T(h)) und (¹₄h², T(¹₂h)), ergibt sich P(T(√

· ) |h², 1

4h²)(x) = T(h) + T(¹₂h) − T(h)

1

4h² − h² (x − h²).

Da man T(0) ann¨ahern will, extrapoliert man an der Stelle x = 0, d.h.

- 6

r

r

0 ¹

4h² h² x

T₁(h) T(¹₂h) T(h)

T(√ x)

Beispiel 10.9.

Sei

I =

Z _π/2

0 xcos x + e^x dx = π

2 + e^12π − 2

und T(h) die zugeh¨orige Trapezregel. Die Extrapolation angewandt auf die Trapezregel liefert folgende Resultate:

n T(h) T₁(h) = ⁴₃T(h) − ¹₃T(2h) |T₁(h) − I| 4 4.39692773

8 4.38523920 4.38134302 6.93e−05

16 4.38226833 4.38127803 4.33e−06

32 4.38152257 4.38127398 2.70e−07

Die N¨aherung T₂(h) l¨aßt sich wie T₁(h) ebenfalls ¨uber Extrapolation erkl¨aren. Es gilt

P(T(√

· ) |h², 1

4h², 1

16h²)(0) = T₂(h).

Allgemeine Vorgehensweise. Sei

T_i,0 := T(2⁻ⁱh), i = 0,1,2, . . . , wobei h eine feste Anfangsschrittweite ist.

Es soll das Interpolationspolynom P(T(√

· ) |h², . . . ,(2^−kh)²)(x) an der Stelle x = 0 ausgewertet werden.

Anwendung des Neville-Aitken Schemas um den Wert T_k(h) = T_k,k = P(T(√

· ) |h², . . . ,(2^−kh)²)(0) zu berechnen, liefert die Rekursion

T_i,j = 4^jT_i,j−1 − T_i−1,j−1

4^j − 1 , j = 1,2, . . . , i ≥ j.

Romberg-Schema

T(h) = T_0,0

ց

T(¹₂h) = T_1,0 → T_1,1

ց ց

T(¹₄h) = T_2,0 → T_2,1 → T_2,2

ց ց ց

T(¹₈h) = T_3,0 → T_3,1 → T_3,2 → T_3,3

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

T_i−1,j−1

ց

−1 4^j−1

T_i,j−1 −→ T_i,j

4^j 4^j−1

Beispiel 10.10.

Sei

I =

Z _π/2

0 xcos x + e^x dx = π

2 + e^12π − 2,

wie in Beispiel 10.2, und T_i,0 = T(2⁻ⁱh), wobei T(·) die Trapezregel ist. F¨ur die Anfangsschrittweite h = ¹₄^π₂ ergibt das Romberg-Schema folgende Werte:

i T_i,0 T_i,1 T_i,2 T_i,3

0 4.396927734684

1 4.385239200472 4.381343022401

2 4.382268326301 4.381278034910 4.381273702411

3 4.381522565173 4.381273978130 4.381273706768 4.381273707762

Fehler:

i |I − T_i,j|

0 1.57e−02

1 3.97e−03 6.93e−05

2 9.95e−04 4.33e−06 5.35e−09

3 2.49e−04 2.70e−07 8.22e−11 1.42e−12

10.5 Zweidimensionale Integrale 10.5.1 Transformation von Integralen.

Eindimensionales Integral

Z _b

a f(x) dx.

Sei

I₁ = [a, b], I₂ = [c, d] und ψ : I₁ → I₂ eine stetig differenzierbare bijektive Abbildung.

Es gilt die Transformationsformel

Z

I₁ f(ψ(x)) ψ^′(x) dx =

Z

I₂ f(y) dy.

Ein interessanter Spezialfall ergibt sich, falls ψ affin ist, d.h.

Wenn Q_m(g;I₁) = (b − a)

m X i=0

w_ig(x_i)

eine Formel zur Ann¨aherung von ^R_a^b g(x) dx ist, ergibt sich eine ent- sprechende Quadraturformel f¨ur das Intervall I₂ = [c, d]:

Z

I₂ f(y) dy =

Z _b

a f( ˆψ(x))|ψˆ^′(x)|dx

= d − c b − a

Z _b

a f( ˆψ(x)) dx ≈ (d − c)

m X i=0

w_if( ˆψ(x_i)) , also insgesamt:

Q_m(g; I₁) = (b − a)

m X i=0

w_ig(x_i) Q_m(f; I₂) = (d − c)

m X i=0

wˆ_if(ˆx_i), mit wˆ_i = w_i, ˆx_i = x_i − a

b − a d + b − x_i b − a c .

Beispiel 10.11.

Gauß-Quadraturformeln werden oft f¨ur das Intervall [−1, 1] spezifiziert, z.B.

Z ₁

−1 f(x)dx ≈ 2

1

2f − 1 3

√3 + 1

2f1 3

√3

aus Beispiel 10.8.

Die entsprechende Formel f¨ur ein beliebiges Intervall [c, d]:

Z _d

c f(x) dx ≈ h 2

fc + (1 2 −

√3

6 )h + fc + (1 2 +

√3 6 )h

, h := d − c.

Analog kann man f¨ur die Gauß-Quadratur mit 4 St¨utzstellen

Z ₁

−1 f(x) dx ≈ 2

3 X i=0

w_if(x_i),

w₀ = w₃ = 0.173928, w₁ = w₂ = 1

2 − w₀,

Wir betrachten nun die Transformation eines zweidimensionalen Integrals

Z

B f(x, y) dx dy, B ⊂ R²_.

Sei B₁, B₂ ⊂ R² _{und ψ : B1 → B2} eine stetig differenzierbare bijektive Abbildung mit Jacobi-Matrix

J(x, y) =





∂ψ₁

∂x (x, y) ^∂ψ_∂y¹(x, y)

∂ψ₂

∂x (x, y) ^∂ψ2

∂y (x, y)



. Es gilt folgende Verallgemeinerung von (10.38):

Satz 10.12. Falls det J(x, y) 6= 0 f¨ur alle (x, y) ∈ B₁, so gilt

Z

B₁ f(ψ(x, y))|detJ(x, y)| dx dy =

Z

B₂ f(˜x,y)˜ d˜x d˜y.

F¨ur den Spezialfall, daß ψ affin ist, ψ(x, y) = Ax

y

+ b, A ∈ R^2×2_{, det(A) 6= 0, b ∈} R²_,

ergibt sich daraus die Transformationsformel

|detA|

Z

B₁ f(Ax y

+ b)dx dy =

Z

B₂ f(˜x,y˜) dx d˜ y.˜

Mit Hilfe dieser Transformationsformel kann man, wie im eindimensio- nalen Fall, eine Quadraturformel f¨ur einen Standardbereich (z.B. Einheitsquadrat, Einheitsdreieck) in eine Formel f¨ur einen

Wichtiger Unterschied zwischen ein- und mehrdimensionaler Integra- tion:

Zwei Intervalle [a, b] und [c, d] lassen sich stets durch affine Transforma- tionen aufeinander abbilden. Hingegen ist es meistens nicht m¨oglich, einfache Gebiete in Rⁿ _{, n ≥} 2, durch eine affine Transformation inein- ander zu ¨uberf¨uhren.

Beispiel 10.14.

Sei B₁ = [0,1] × [0,1] das Einheitsquadrat.

Jede affine Abbildung bildet B₁ auf ein Parallelogramm ab. Eine affine Abbildung von B₁ auf den Einheitskreis S = {(x, y)|(x² + y²) ≤ 1} ist

also nicht m¨oglich. △

Affine Transformationen

- 6

B₁

p

p p

p

B₂

0 1 2 3 4 5 6

1 2 3 4

_ψ 3

- 6

@@

B₁ @

p

p p

((((XX((((XXX(X

B₂

0 1 2 3 4 5 6

1 2 3

1

ψ

10.5.2 Integration ¨uber dem Einheitsquadrat

Z ₁ 0

Z ₁

0 f(x, y) dx dy.

Sei

Q_m(g) =

m X i=0

w_ig(x_i)

eine Quadraturformel f¨ur das eindimensionale Integral ^R₀¹ g(x) dx und F(y) :=

Z ₁

0 f(x, y) dx

Z ₁ 0

Z ₁

0 f(x, y) dx dy =

Z ₁

0 F(y) dy ≈

m X j=0

w_jF(x_j)

=

m X j=0

w_j

Z ₁

0 f(x, x_j) dx ≈

m X j=0

w_j

m X i=0

w_if(x_i, x_j)

=

m X i,j=0

w_iw_j f(x_i, x_j) =: Q⁽²⁾_m (f).

Beispiel 10.17.

Sei

Q₁(g) = 1

2g(x₀) + 1

2g(x₁), x₀ := 1 2 −

√3

6 , x₁ := 1 2 +

√3 6 ,

die eindimensionale Gauß-Quadraturformel mit zwei St¨utzstellen f¨ur das Intervall [0,1].

Daraus ergibt sich die Produktregel Q⁽²⁾₁ (f) = 1

4f (x₀, x₀) + 1

4f (x₀, x₁) + 1

4f (x₁, x₀) + 1

4f (x₁, x₁)

f¨ur den Bereich [0,1] × [0,1]. Diese Formel ist exakt f¨ur alle Linear-

10.5.3 Integration ¨uber dem Einheitsdreieck F¨ur Dreiecke ist es zweckm¨aßig, von den Monomen

1, x, y, x², xy, y², usw.

auszugehen und die Frage nach solchen Quadraturformeln zu stellen, die alle Monome der Form x^k1y^k2, 0 ≤ k₁ + k₂ ≤ M exakt integrieren.

Einige typische Beispiele:

(i) Q(f) = 1

2f(1 3, 1

3) (ii) Q(f) = 1

6[f(0,0) + f(1,0) + f(0,1)]

(iii) Q(f) = 1

6[f(1

2,0) + f(0, 1

2) + f(1 2, 1

2)]

(iv) Q(f) = 1

6[f(1 6, 1

6) + f(2 3, 1

6) + f(1 6, 2

3)].

Die Monome 1, x, y werden durch die Formeln in (i), (ii) exakt inte- griert.

Die Monome 1, x, y, xy, x², y² werden durch die Formeln in (iii), (iv) exakt integriert.

Kapitel 11 Gew¨ohnliche Differentialgleichungen

11.1 Einf¨uhrung

Gesucht wird eine Funktion y = y(t) einer (Zeit-)Variablen t, die der Gleichung

y^′(t) = f(t, y(t)) , t ∈ [t₀, T] , und der Anfangsbedingung

y(t₀) = y⁰ gen¨ugen soll.

Man spricht von einer gew¨ohnlichen Differentialgleichung erster Ord-

Allgemeiner hat man mit Systemen von n gew¨ohnlichen Differential- gleichungen erster Ordnung

y₁^′ (t) = f₁(t, y₁(t), . . . , y_n(t)) ...

y_n^′ (t) = f_n(t, y₁(t), . . . , y_n(t)) zu tun. Anfangsbedingungen:

y_i(t₀) = y_i⁰, i = 1, . . . , n, Setzt man

y(t) :=







y₁(t) ...

y_n(t)





, f(t, y) :=







f₁(t, y₁(t), . . . , y_n(t)) ...

f_n(t, y₁(t), . . . , y_n(t))





, y⁰ :=







y₁⁰ ...

y_n⁰





, kann dieses Anfangswertproblem kompakt wieder in der Form

y^′(t) = f(t, y(t)) , t ∈ [t₀, T] , y(t₀) = y⁰ geschrieben werden.

Beispiel 11.2.

Gesucht werden Funktionen y₁(t), y₂(t), t ≥ 0, f¨ur die

y₁^′ y₂^′

!

=

1

2y₁ − y₂

2y₁ − 2y₂ + 3 sint

!

(t ≥ 0) und y₁(0) y₂(0)

!

= 2 1

!

gilt. Durch Einsetzen kann man einfach nachpr¨ufen, daß y₁(t)

y₂(t)

!

= 2 cos t cos t + 2 sint

!

die L¨osung ist. △

Beispiel 11.4.

Chemische Reaktionsprozesse werden oft mit gew¨ohnlichen Differen- tialgleichungen modelliert. Seien S₁, . . . , S_n chemische Stoffe, die bei konstanter Temperatur in einem abgeschlossenen System mit einander reagieren. Die i-te Reaktion wird durch

n X j=1

a_ijS_j −→^ki

n X j=1

b_ijS_j

beschrieben. Als Beispiel betrachten wir die chemische Pyrolyse S₁ −→^k1 S₂ + S₃

S₂ + S₃ −→^k2 S₅ S₁ + S₃ −→^k3 S₄

S₄ −→^k4 S₃ + S₆.

Sei y_i(t) die Konzentration der Komponente S_i zum Zeitpunkt t. Das zugeh¨orige System gew¨ohnlicher Differentialgleichungen, das die Dy- namik dieses Reaktionsprozesses beschreibt, ist

y₁^′ = −k₁y₁ − k₃y₁y₃ y₂^′ = k₁y₁ − k₂y₂y₃

y₃^′ = k₁y₁ − k₂y₂y₃ − k₃y₁y₃ + k₄y₄ y₄^′ = k₃y₁y₃ − k₄y₄

y₅^′ = k₂y₂y₃ y₆^′ = k₄y₄.

△

Beispiel 11.5.

Die Entwicklung der Temperatur T(x, t) an der Stelle x eines Stabes zur Zeit t ergibt sich als L¨osung der Anfangsrandwertaufgabe

∂T

∂t = κ∂²T

∂x² , t > 0, x ∈ (0, ℓ).

Die Anfangswerte seien T(x,0) = Φ(x). Die Randwerte sind T(0, t) = T(ℓ, t) = 0.

Mit Hilfe der sogenannten Linien-Methode kann man die gesuchte L¨osung T(x, t) mit Hilfe eines Systems gew¨ohnlicher Differentialglei- chungen erster Ordnung ann¨ahern.

κ∂²T(x, t)

∂x² ≈ κ

h²_x (T(x + h_x, t) − 2T(x, t) + T(x − h_x, t)), y_j(t) ≈ T(x_j, t), j = 1,2, . . . , n_x − 1,

y_j^′ = κ

h²_x(y_j+1 − 2y_j + y_j−1), j = 1, . . . , n_x − 1,

Mit y = (y₁, y₂, . . . , y_nx−1)^T ergibt sich

y^′ = f(t, y) = Ay, wobei A die Tridiagonalmatrix

A = − κ h²_x

















2 −1

−1 2 −1 ∅ . . . .

∅ −1 2 −1

−1 2

















ist. △

Gew¨ohnliche Differentialgleichung m-ter Ordnung: in der Differential- gleichung kommen Ableitungen bis zur Ordnung m vor.

Beispiel 11.6. Das mathematische Pendel in Beispiel 1.2 wird durch die gew¨ohnliche Differentialgleichung zweiter Ordnung

φ^′′(t) = −g

ℓ sin(φ(t)) und die Anfangsbedingungen

φ(0) = φ₀, φ^′(0) = 0

charakterisiert. △

Anfangswertaufgabe m-ter Ordnung: Bestimme eine skalare Funktion y(t), so daß

y^(m) = g(t, y, y^′, . . . , y^(m−1)) , t ∈ [t₀, T] ,

y(t₀) = z₀, y^′(t₀) = z₁, . . . , y^(m−1)(t₀) = z_m−1.

11.2 Reduktion auf ein System 1. Ordnung Setzt man

y₁(t) := y(t)

y₂(t) := y^′(t) = y₁^′ (t) y₃(t) := y^′′(t) = y₂^′ (t)

...

y_m(t) := y^(m−1)(t) = y_m^′ ₋₁(t), so folgt aus (11.10) y_m^′ = g(t, y₁, y₂, . . . , y_m), also

y₁^′ (t) = y₂(t) y₂^′ (t) = y₃(t)

...

y_m−^′ ₁(t) = y_m(t)

y_m^′ (t) = g(t, y₁(t), . . . , y_m(t))



















f¨ur t ∈ [t₀, T]

Beispiel 11.8.

Die gew¨ohnliche Differentialgleichung dritter Ordnung y^′′′ = −2y^′′ + y^′ + y² − e^t, t ∈ [0, T], mit Anfangsbedingungen

y(0) = 1, y^′(0) = 0, y^′′(0) = 0, ergibt ¨uber

y₁(t) := y(t), y₂(t) := y^′(t), y₃(t) := y^′′(t), das ¨aquivalente System erster Ordnung







y₁^′ y₂^′ y₃^′





 =







y₂ y₃

−2y₃ + y₂ + y₁² − e^t





 (t ∈ [0, T]) mit Anfangsbedingungen

(y₁(0), y₂(0), y₃(0)) = (1, 0,0). △

11.3 Einige theoretische Grundlagen

Ein Problem heißt korrekt gestellt, falls

1.) eine L¨osung existiert,

2.) diese L¨osung eindeutig ist und

3.) stetig von den Daten abh¨angt.

Nach dem Satz von Peano existiert bereits unter sehr schwachen An- forderungen an die rechte Seite f, n¨amlich Stetigkeit in t und y, eine

Unter etwas st¨arkeren Voraussetzungen an f gilt allerdings der fol- gende Existenz- und Eindeutigkeitssatz von Picard-Lindel¨of , der die Forderung 2.) sichert.

Satz 11.10. Sei f eine Funktion, die stetig in (t, y) und dar¨uber hinaus im folgenden Sinne Lipschitz-stetig in y ist (wobei k · k eine beliebige feste Norm auf Rⁿ _ist):

kf(t, y) − f(t, z)k ≤ Lky − zk

f¨ur alle t ∈ [t₀, t₀ + δ] mit δ > 0 und f¨ur alle y, z aus einer Umge- bung U von y⁰.

Dann existiert eine eindeutige L¨osung y von (11.3) in einer Um- gebung von t₀ (die von δ, kfk und von U abh¨angt).

Beispiel 11.11.

Wir w¨ahlen k · k = k · k_∞, die Maximumnorm. F¨ur das mathematische Pendel aus Beispiel 11.7 ergibt sich mit c := −g/ℓ

kf(t, y) − f(t, z)k_∞ =

y₂ csiny₁

!

− z₂ csin z₁

! _∞

=

y₂ − z₂

ccos(ξ)(y₁ − z₁)

! _∞

= max{|y₂ − z₂|,|c| |cos(ξ)| |y₁ − z₁|}

≤ max{1,|c|}ky − zk_∞,

also gilt die Lipschitz-Bedingung mit L = max{1, ^g} f¨ur alle t, y, z. △

Forderung 3.) der Korrektgestelltheit ist nat¨urlich f¨ur numerische Zwecke besonders essentiell.

Im vorliegenden Zusammenhang betrachten wir den einfachsten Fall, daß unter

”Daten“ lediglich die Anfangswerte y⁰ verstanden werden.

Satz 11.12. Die Funktion f erf¨ulle die Lipschitz-Bedingung (11.12) (bzgl. einer Umgebung U von y⁰, z⁰ ∈ Rⁿ_).

Seien y(t), z(t) L¨osungen von (11.3) bez¨uglich der Anfangsdaten y⁰, z⁰ ∈ Rⁿ_.

Dann gilt f¨ur alle t aus einer Umgebung von t₀ die Absch¨atzung ky(t) − z(t)k ≤ e^L|t−t0|ky⁰ − z⁰k.

Bemerkung 11.14.

Die Funktion y(t) l¨ost

y^′(t) = f(t, y(t)), y(t₀) = y⁰, genau dann, wenn sie

y(t) = y⁰ +

Z _t

t₀ f(s, y(s))ds l¨ost.

11.4 Einfache Einschrittverfahren

Man verwendet y¹ = y⁰ + hf(t₀, y⁰) als N¨aherung f¨ur y(t₀ + h).

- 6

y⁰

t₀ t₀ + h

y(t)

y¹ q q

Dies f¨uhrt auf das folgende Verfahren:

Algorithmus 11.16 (Euler-Verfahren) Gegeben: Schrittweite h = ^T−t0

n mit n ∈ N. Berechne f¨ur j = 0, . . . , n − 1:

t_j+1 = t_j + h

y^j+1 = y^j + hf(t_j, y^j).

Ein Herleitungsprinzip

Bemerkung 11.17. Sei die Anfangsbedingung (t_j, y^j) ∈ R² _gegeben.

Die Funktion ˜y(t) l¨ost das lokale Anfangswertproblem y˜^′ = f(t,y)˜ f¨ur t ∈ [t_j, T], y(t˜ _j) = y^j, genau dann, wenn

y(t) =˜ y^j +

Z _t

t_j f(s,y(s))˜ ds, t ∈ [t_j, T], gilt. Insbesondere gilt f¨ur t = t_j+1:

y(t˜ _j+1) = y^j +

Z _tj+1

t_j f(s,y(s))˜ ds.

△

Das Euler-Verfahren erh¨alt man dann ¨uber