Gewöhnliche Differentialgleichungen - tu-ilmenau.de · Viele praktische Beispiele führen auf...

Technische Universitat Ilmenau Postfach 10 0565

Fakultat fur Mathematik 98684 Ilmenau

und Naturwissenschaften Germany

Institut fur Mathematik Tel.: 03677/692652

Fax: 03677/691241

Telex: 33 84 23 tuil d.

email: [email protected]

Preprint No. M 11/98

Gewohnliche

Differentialgleichungen

Beispiele, Modelle, Verfahren, Software

Werner Neundorf

Mai 1998

‡MSC (1991): 58F40, 65-01, 65-04, 65L05, 68Q25

Viele praktische Beispiele fuhren auf die Anwendung von Differentialgleichungen. Ausder Literatur wurden zahlreiche Modelle zitiert und erlautert. Ihre Analyse fallt indas Gebiet der Theorie der gewohnlichen Differentialgleichungen und findet hier kei-nen Niederschlag. Zur Numerik der GDGL wird bis auf wenige Verfahrensvorschlageauf die einschlagige Literatur verwiesen. Zusatzlich findet man einige Hinweise aufentsprechende Software.

Many practical examples lead to the application of differential equations. Startingfrom the literature, it goes on to present and describe a lot of models. Their mosttheoretical analysis comes from the Theory of Ordinary Differential Equations andwas not involved here. The techniques for the numerical solution of ODE are presentin some problems, but the reader finds a large presentation of numerics in manytextbooks. In addition, there are given some software hints.

Key words: applications for ordinary differential equations, tutorial aspects,analysis of algorithms, programs

MSC (1991): 58F40, 65-01, 65-04, 65L05, 68Q25

1 Einleitung

Teile dieser Arbeit basieren auf einem Skriptum, daß ich in der LehrveranstaltungNumerik der gewohnlichen Differentialgleichungen fur die Studenten der Mathematikim 5. Semester sowie der Allgemeinen und Theoretischen Elektrotechnik einbezogenhabe. Dabei sollte ein Uberblick zur Darstellung und Simulation verschiedenartigerSachverhalte durch Differentialgleichungen gegeben werden. Entsprechende prakti-sche Anwendungen auf den Gebieten Technik, Physik, Chemie, Mathematik, Biolo-gie, Okonomie, Okologie und Medizin sind aus meistens einschlagiger Literatur sowieSammlungen von Angehorigen des Fachgebiets Numerische Mathematik und Infor-mationsverarbeitung zusammengetragen worden.Die Aufgabenstellungen bieten auch eine breite Auswahl fur mogliche Beleg- undPraktikumsarbeiten mit programmtechnischer Umsetzung.Bezuglich der Theorie der gewohnlichen Differentialgleichungen mochte ich insbeson-dere auf das Skriptum von J.Knobloch und J.Steigenberger zur Vorlesung Gewohnli-

che Differentialgleichungen fur Studenten der Mathematik im 4. Semester im Rahmender Grundausbildung in Analysis verweisen.

2

2 Differentialgleichungen in der Praxis

2.1 Elektrotechnik

1. Elektrische Schaltung [Schwetlick, Kretzschmar]

Aus der Beziehung fur den Spannungsabfall berechnet man die Stromstarke.

cc

r r∼

R L

IU

U(t) = U0 sin(ω0t), t0 = 0 ≤ t ≤ T

I(t0) = 0

DGL: U(t) = R I(t)︸︷︷︸

Ohmscher

+ LdI(t)

dt︸︷︷︸

Induktiver

Spannungsabfall

Allgemeine Losung fur die Stromstarke:

I(t) = e−∫

R

Ldt

[

c0 +

∫R

Le∫

R

Ldτdt

]

, c0 aus Anfangsbedingung ermitteln.

Losung im speziellen Fall:

R = R0 = const, L = L0 = const, U(t) = U0 sin(ω0t), U0 = const, t0 = 0

I(t) =U0ω0L0

R20 + ω2

0L20

e−

R0

L0t+

U0√

R20 + ω2

0L20

sin(ω0t− α) mit α = arctanω0L0

R0

2. Einschalten eines Stromkreises [Philippow]

Anwendung der Kirchhoffschen Gesetze.

DGL:dU

dt+

1

RCU =

1

RCU0, I(t) = C

dU(t)

dt

Transformation der Variablen liefert eine DGL 1.Ordnung

x+ f x = g(τ), τ transformierte Zeit

3

-

r rUo

R C

I

U

Hangen C und/oder R auch vom Strom I ab (nichtlineare Kapazitaten oderWiderstande), so gilt

x+ f(x) x = g(x, τ).

Diese DGL ist i.a. nicht geschlossen integrierbar.

3. Nichtlinearer ungedampfter Schwingkreis [Philippow]

Man betrachtet man den ungedampften Schwingkreis mit nichtlinearer Induk-tivitat.

r rC L

I

Wenn der Zusammenhang zwischen Strom I(t) und dem Fluß ψ(t) die KennlinieI = aψ + bψ3 besitzt, so erhalt man die DGL

d2ψ

dt2+a

Cψ +

b

Cψ3 = 0.

Nach Normierung mit Ψ = maxψ(t) = ψ(0) durch

x(t) =ψ(t)

Ψ, τ =

√a

Ct, λ = Ψ2a

b

erhalt man die DGL 2.Ordnung

x+ x+ λx3 = 0,

die nicht geschlossen losbar ist.Das Anlegen einer sinusformigen Erregerschwingung an diesen Schwingkreis lie-fert eine erzwungene Schwingung mit der normierten DGL (Duffingsche DGL)

x+ x+ λx3 = k cos(ατ) .

4

4. Relaxationsschwingungen [Philippow]

Zahlreiche Probleme der nichtlinearen Elektrotechnik (ungedampfte Schwin-gungen, Vibrator, Multivibrator) fuhren zu der normierten DGL 2.Ordnung

y − ε (1− y2) y + y = 0, t ≥ 0, ε ≥ 0,

(Van der Pol DGL).

Fur ε = 0 liefert sie bei Anfangsbedingungen y0, y0 an der Stelle t = 0 dieharmonische Schwingung y(t) = y0 sin t+ y0 cos t mit der Periode T = 2π.Wachst ε>0 an, so verzerrt sich die Schwingung immer mehr und die Periodevergroßert sich.

Bei ε≫ 1 (z.B ε = 100, 1000) beschreibt die DGL Relaxationsschwingungen.Dabei wechseln sich Zeiten sehr langsamer Bewegung mit solchen sehr schnellerBewegung ab.Der Parameter ε legt auch die Einschwingdauer fest.Bei erregten Schwingkreisen steht ein zusatzlicher Term auf der rechten Seiteder DGL.

5. Mehrmaschige elektrische Netzwerke [Philippow, Uhlmann]

Das zweimaschige gekoppelte Netzwerk fur die normierten Variablen x1(t), x1(t)hat die Form

x1 + α1x1 + ε1(1− x22)x2 + x1 = 0x2 + α2x2 + ε2(1− x21)x1 + x2 = 0

Eine geschlossene Losung ist nicht moglich. Insbesondere sind Losungen furgroße Werte von |ε1|, |ε2| gesucht.

6. Dynamische Kennlinie des Asynchronmotors [Goering]

Bei Asynchronmotoren ist die Drehzahl stark von der Belastung abhangig. Furschwankende Motorleistung ist es zur Berechnung von Torsionsschwingungennotwendig, den Zusammenhang zwischen Motormoment und Winkelgeschwin-digkeit zu bestimmen.Mit normierten Großen u(t), v(t), t dimensionslose Zeit, erhalt man das nichtgeschlossen losbare DGL-System 1.Ordnung

βu = x− u+ xv, u(0) = u0, β ∈ R+, x = x(t)

βv = −xu − v, v(0) = v0.

7. Weitere Beispiele

(a) Spannungsstabilisatorschaltung [Philippow]

Ferroresonanzstabilisator, y = y(t)

y − δy + ay + by9 + αy = −Γ cos(t+ φ)

5

(b) Frequenzteilerschaltung

y − ε(1− y2 − y2)y + y + 2Byy = B sin(2t)

(c) System der Drei-Wellen-Wechselwirkung

a0 = a0 + 2a0a21 sin ϑ

a1 = −a1 − a20a1 sin ϑ

a2 = 2(a21 − a20)− 2δ + 2(2a21 − a20) cos ϑ

(d) Rossler-System

X = − Y − Z

Y = X + aY

Z = bX − cZ +XZ

2.2 Mechanik

1. Vibration von Membranen und Staben [Holmes, Seydel]

Diese fuhren auf DGL vom Duffingschen Typ

y − ε y − βy + αy2 = f cos(ωt)

Sie stellt eine erste Naherung fur die Schwingung dar und wird fur verschiedeneParameterkonstellationen ε, β, α, f, ω analysiert.

2. Schwingungen von Eisenbahnwagen [Kubicek, Marek]

Ix+ 2lfT

(lx

v+r1 − r22r0

)

+Kxx = 0

My + 2fL

(y

v− x

)

Kyy + Fg(y − δ) = 0

3. Bewegungen eines Himmelskorpers [Shampine, Gordon]

Die Bewegung einer Rakete oder eines Satelliten im Schwerefeld zweier großerHimmelskorper fuhrt auf ein restringiertes Dreikorperproblem.Bsp.: Erde - Mond - Rakete

Sonne - Erde - RaketeSonne - Saturn - Saturnmond

y′′1 = 2y′2 + y1 −µ∗(y1 + µ)

r31− µ(y1 − µ∗)

r32

y′′2 = −2y′1 + y2 −µ∗y2r31

− µy2r32

mit µ∗ = 1−µ, 0 < µ < 1 und r1 =√

(y1 + µ)2 + y22, r2 =√

(y1 − µ)2 + y22.

6

4. Orientierung von Telekommunikationssatelliten [Troger]

Dabei geht es um die positive Orientierung der Antenne des Satelliten auf seinerelliptischen Umlaufbahn durch die Vorgabe einer geeigneten Anfangsorientie-rung.

-

6

u cz

* Yν ϕ

(1 + e cos ν)d2ϕ

dν2− 2e sin ν

dϕ

dν+ α sinϕ cosϕ = 2e sin ν

e Extentrizitat, α Konstante zur Beschreibung der Massenverteilung,ν Umlaufwinkel, ϕ gesuchte Abweichung vom Leitstrahl

5. Pendel mit Reibung [Arnold]

x+ kx+ x = 0, x(0), x(0) gegeben, k ≥ 0 Reibungskoeffizient.

Uberfuhrung in ein DGL-Systemx = y, x(0) gegeben

y = x = −kx− x = −x− ky, y(0) = x(0) gegeben

Notation in Matrixform

X = AX =

(0 1-1 -k

)(xy

)

, trace(A) = −k, det(A) = 1

Eigenwertproblem zur Matrix A

λ(A) : 0 = det(A− λI) = λ2 + kλ+ 1.Falls |k| > 2, dann λi reell, negativ und verschieden.

Sei

Λ =

(λ1 00 λ2

)

Eigenvektoren

x1 =

(-31

)

, x2 =

(1-3

)

, U = (x1, x2) regular, AU = UΛ

7

Transformation auf Diagonalform moglich.

X = UZ

X = UZ = AX = AUZ

UZ = UΛZ

Z = ΛZ

Z1 = λ1Z1, entkoppeltes System

Z2 = λ2Z2

Z1(t) = eλ1tZ1(0), Losungen

Z2(t) = eλ2tZ2(0)

Beispiel

X0 = (x(0), y(0))T = (2, 2)T , k = 103, λ1 = −1

3, λ2 = −3

1. Darstellung der Losungen als Phasenkurven im Raum der EV

-

6

Z1

Z2

?Rj

3

6 I

k

λ2 < λ1 < 0

2. Darstellung der Losungen als Phasenkurven im Raum (x1, x2)

U−1 = 18

(-3 -1-1 -3

)

, Z0 = U−1X0 =

(-1-1

)

X = UZ =

(-3 11 -3

)(eλ1tZ1(0)eλ2tZ2(0)

)

=

(3e−t/3 − e−3t

−e−t/3 + 3e−3t

)

Bei gegebener Anfangsauslenkung kommt es durch den Anstoß zu einem kurz-zeitigen Uberschwingen des Pendels. Danach wird die Pendelbewegung wegender Reibung jedoch stark gedampft.Die Losung x(t) = 3e−t/3 − e−3t strebt mit wachsender Zeit asymptotischgegen Null.

8

-

6

x

y

x1

x2

s

s

w

o

K

U

K

U

2.3 Thermodynamik, Chemie

1. Lorenzmodell der Warmestromung in einer Flussigkeits-schicht [Lorenz, Kubicek, Steeb]

Eine Flussigkeitszelle wird von einer Seite erwarmt und von der gegenuber-liegenden Seite gekuhlt. Die dabei resultierende konvektive Bewegung ist vonInteresse.Das DGL-System entsteht aus den Navier-Stokes-Gleichungen.

y1 konvektive Durchmischungy2 horizontale Temperaturveranderungy3 vertikale Temperaturveranderung

y′1 = −σy1 + σy2, y1(0) = y10

y′2 = −y1y3 + ry1 − y2, y2(0) = y20

y′3 = y1y2 − by3, y3(0) = y30

σ Prandtl-Zahl, r reduzierte Rayleigh-Zahl, b Wellenzahl sind positive Para-meter.

2. Schwingungen im Warmeaustauscher [Troger]

Durch Warmeunterschiede in einem solchen Durchlaufsystem konnen Rohrenunter Spannung, andere unter Druck stehen. Dadurch entstehen auch Schwin-gungen.

y + ε1y + ε2y + y3 + y3 = 0

9

3. Chemische Reaktoren [Kubicek, Marek]

Beschreibung einer Kaskade von 2 chemischen Tankreaktoren mit Ruhrwerkund Ruckkopplung liefert ein System von 4 DGL 1.Ordnung.

- - -∞ ∞y1T1 y2T2

y1

T1

y2

T2

1− c

yi Konzentration, Ti Temperatur, 0 < c ≤ 1

2.4 Physik

1. Freier Fall

Angenommen, das AWPy′ = a

√y, y = y(t), y(0) = c, a = const > 0, c = const ≥ 0,

beschreibt den freien Fall eines Korpers, y(t) sei also der zur Zeit t zuruckgelegteWeg.

(a) Welche Losung erhalt man im Fall c = 0?Gibt es noch eine zweite Losung der DGL mit der AB y(0) = 0?

(b) Mit welcher AB kann ein Einschrittverfahren einen freien Fall berechnen?

(c) Gegeben sei τ ≥ 0 und die Funktion

z(t) =

0 , falls 0 ≤ t ≤ τ

14c2(t− τ)2 , falls t > τ.

Erfullt z(t) die DGL?Welche Auswirkung hat die Antwort auf das Problem der Eindeutigkeitder Losung?Worin liegt die Ursache fur diese Situation?

Hinweis: Untersuche die Glattheit der rechten Seite f(t, y).

2. Bewegungsgleichung eines Korpers in laminarer Stromung

Die Bewegungsgleichung eines Korpers mit der Masse m in einer zahen Flussig-keit (laminare Stromung) mit dem Stokes’schen Widerstandgesetzes R(v) =K1v (K1 Materialeigenschaften des Mediums/Geometrie des Korpers) ergibtsich aus dem Kraftgesetz

10

mv = A+G+R(v), G Erdanziehung, A Auftrieb.

Die abgeleitete Modellgleichung fur die Geschwindigkeit v(t) sei

v = c + k1v, c > 0, v(0) gegeben.

- Diskutiere eine Losung bei verschiedenen Startwerten v(0).

- Es sei v(0) = 0, k1 < 0, c > 0.Gibt es eine Gleichgewichtsgeschwindigkeit v∞ = v(t = ∞) = const?

3. Bewegungsgleichung eines Korpers in turbulenter Stromung

Die Bewegungsgleichung eines Korpers mit der Masse m in einer turbulentenStromung (Neigung zur Wirbelbildung) mit dem Newtonschen Widerstandge-setzes R(v) = K2v

2 (K2 Materialeigenschaften/Geometrie) ergibt sich aus demKraftgesetz

mv = A+G+R(v).

Die abgeleitete Modellgleichung fur die Geschwindigkeit v(t) seiv = c + k2v

2, c > 0.

Es sei v(0) = v0, k2 < 0, c > 0.- Diskutiere eine Losung bei verschiedenen Startwerten.- Gibt es eine Gleichgewichtsgeschwindigkeit v∞ = v(t = ∞) = const?

4. Bewegungsgleichung fur den freien Fall

Die Bewegungsgleichung fur den freien Fall (ohne Reibung) mit der von derHohe z(t) abhangigen Erdbeschleunigung g(z) lautet z = −g(z).Nach dem Newton’schen Gravitationsgesetz gilt g(z) =

g0z20

z2, wobei g0 = g(z0)

den Wert der Erdbeschleunigung auf der Hohe z0 bezeichnet.Fur kleine Fallhohen kann man 3 Approximationen der rechten Seite der DGLanwenden.

(1) g0 = g = 9.806 ms−2 Normwert, rE = 6.375 · 106 m Erdradius, z0 = z,

z = −g,(2) g(z) ≈ g − 2

g

rE(z − rE) = 3g − 2

g

rEz,

z = −3g + 2g

rEz,

(3) z0 = rE, g0 = g,

z = −g r2E

z2.

Die AB sind jeweils z(0) = rE +H, z(0) = 0, H > 0 Fallhohe.

(a) Welche konstante Funktion ist Losung der Modellgleichung (2)?

11

Finde ihre allgemeine Losung.

(b) Berechne die Fallzeit tF und Endgeschwindigkeit z(tF ) bei einer FallhoheH und gegebenen AB fur die Modellgleichungen (1), (2).

(c) Welchen Einfluß hat die Linearisierung des Modells?

5. Bewegungsgleichungen, Schwingungen, Pendel(parametrischen DGL)

(a) Schwingungsgleichung ohne Einwirkung außerer Krafte

y′′ + 2hy′ + k2y = 0

(b) DGL der erzwungenen Schwingung y(t) mit Einwirkung außerer KrafteUntersuche y(t) fur t→ ∞.

y′′ + 2ρy′ + ω2y = a sin(ω1t), y(0) = 0, y′(0) = ω

(c) Van der Polsche DGL

y − ε(1− y2) y + y = 0, y(0) = y0, y(0) = y0, t ≥ 0, ε ≥ 0

Der Parameter ε legt die Einschwingdauer fest.Allgemeine Losung fur ε = 0 : y(t) = y0 sin t+ y0 cos t

Van der Polsche DGL mit periodischer Erregung

y − ε(1− y2) y + y = k cos(Ωt)

(d) Rotation eines Pendels mit der Winkelgeschwindigkeit ω

y + sin y − 1

2ω2 sin(2y) = 0, y(0) =

π

4, y(0) = 0, ω ∈ [0, 10]

(e) Pendel mit Anfangsbedingung

y′′ + k2y = 0, k > 0, y(0) = y0, y′(0) = y′0

Losung: y =y′0ksin(kx) + y0 cos(kx)

y =y′0ksin(kx), falls y0 = 0

y = sin x, falls y0 = 0, y′0 = 1, k = 1

y = c1 sin(kx) + y0 cos(kx), falls AB y′(0) fehlt

(f) Fadenpendel

y′′ + ay′ + by = 0, D = a2 − 4b DiskriminanteLosung:

12

(1) D > 0

y = c1ek1x + c2e

k2x

k1 = −1

2(a−

√D), k2 = −1

2(a+

√D)

z.B. falls b < 0, dann D > 0

(2) D = 0

y = c1ek1x + c2xe

k1x, k1 = −a2

z.B. falls b = 0, a = 0, dann D = 0

(3) D < 0

y = e−ax

2 (c1 sin(12

√

|D|x) + c2 cos(12

√

|D|x)

k1 = −a2+

1

2

√−D i, k2 = −a

2− 1

2

√−D i

z.B. falls b ≫ 0, dann D < 0 ; a > 0 → gedampftes Pendel,a < 0 → erregtes Pendel

(g) Bewegungsgleichung

y+ f(y) = 0, y(0) = 2, y(0) = 0, f(z) =

−1 , falls z ≤ −1z , falls − 1 < z < 11 , falls 1 ≤ z

(h) Duffingsche Bewegungsgleichungen, Vibration y = y(t)

y + y + λy3 = 0, y(0) = y0, y(0) = y0, λ ∈ R+

y + y + λy3 = k cos(αt), λ ∈ R+, k, α ∈ R

y + εy − βy + αy3 = f cos(ωt)

y + 0.7y + y3 = 0.75 cos t, y(0) = y(0) = 0

y + 0.2y + y|y| = 1.5 cos(2t) + 0.5, y(0) = y(0) = 0

y + y − 1

6y3 = 0.8 sin(0.27ωt), y(0) = y(0) = 0, ω = 0.92845

(i) Bewegungsgleichung mit quadratischer Dampfung

y +1

2y2 + y = 0, y(0) = 1, y(0) = 0

(j) Emden’s Bewegungsgleichungend

dt(t2 y) + t2 yn = 0, y(0) = 1, y(0) = 0, n = 2, 3, 5, ...

1

t2d

dt(t2 y) + f(y) = y +

2

ty + f(y) = 0, y(0) = 1, y(0) = 0,

f(y) = sin y bzw. cos y

(k) Liouvillesche DGL

13

d

dt(t2 y) + t2 ey = 0, y(0) = y(0) = 0

(l) Mathieu’s DGL fur ein Pendel mit zeitabhangiger Lange l(t)

y + (a− 2q cos(2t)) y = 0, y(0) = 1, y(0) = 0

a = 1.0, q = 0.1 aufschaukelnde Schwingunga = 1.2, q = 0.1 gedampfte Schwingung

Man kann einen (a, q)-Bereich finden, in dem die Losung beschrankt bleibt(Stutt-Diagramm).

(m) Gedampfte Schwingung einer Membran, Besselsche DGL

t2y′′ + ty′ + (t2 − k2)y = 0, k ∈ Z

Losung : y = J0(t), falls k = 0, y(0) = 1, y′(0) = 0, Besselfunktion1.Art

y = J1(t), falls k = 1, y(0) = 0, y′(0) = 0.5, Besselfunktion1.Art

(n) Resonanz zwischen Eigenschwingung und Storung

y′′ + 6y′ + 5y = 3e−t

(o) Abrollen einer Kugel auf einer Fallinie eines parabolischen Zylinders

y +y

2025 + y2(31532 + y2) = 0, y(0) = 30, y(0) = 0

(p) Rollpendel(Θ

m+ (r2 + s2 − 2rs cos φ)

)

φ+ (rφ2 + g)s sin φ = 0

(q) Spannungsstabilisatorschaltung, Ferroresonanzstabilisator [Philippow]

y − δy + ay + by9 + αy = −Γ cos(t+ φ)

(r) Frequenzteilerschaltung

y − ε(1− y2 − y2)y + y + 2Byy = B sin(2t)

(s) Schwingungen im Warmeaustauscher

y + ε1y + ε2y + y3 + y3 = 0

(t) Bewegungsgleichung fur den freien Fall

z = −gr2E

z2, z(0) = rE +H, z(0) = 0

H Fallhohe, g0 = g = 9.806ms−2 Normwert der Erdbeschleunigung,rE = 6.375 106m Erdradius

14

(u) Hamiltonsche Bewegungsgleichung fur das Henon Heiles-Modell [Steeb]

q1 = p1, p1 = −q1 − 2q1q2

q2 = p2, p2 = −q2 − q21 + q22

AB: z.B. q1(0) = 0.5, q2(0) = p1(0) = p2(0) = 0.1

(v) Bewegungsgleichung fur die biquadratische Hamiltonfunktion

q1 = p1, p1 = −1− r

3q31 − q1q

22, r ∈ [0, 1]

q2 = p2, p2 = −1− r

3q32 − q21q2

Fur wachsendes r wird das System mehr und mehr chaotisch.

(w) Bewegungsgleichung fur die Hamiltonfunktion Toda lattice

aj = aj(bj − bj+1), j = 1, 2 (a3 ≡ 0)

b1 = −2a21b2 = 2(a21 − a22)

b3 = 2a22

(x) Schwingende Atwood-Maschine [Steeb]

(1 + µ) r = r(Θ)2 + cosΘ− µ, µ Masse

r Θ = −2rΘ− sinΘ

(y) Chaotische Systeme mit Orbits (Homoklinen) [Steeb]siehe auch Duffingsche Bewegungsgleichungen

u− 12u+ u2 + u3 = 0

u+ u3 = 0

u+ au+ bu+ cu3 = k1 + k2 cos(Ωt)

Parametersatz z.B. a = 1, b = −10, c = 100, k1 = 0, k2 = 1.2,Ω = 3.5

u+ au+ b sin u = k cos(Ωt)

a = 0.2, b = 1, k = 1.1,Ω = 0.8

u+ au+ (1 + k cos(Ωt)) sin u = 0

a = 0.15, k = 0.94,Ω = 1.56

15

6. Fadenpendel

Wie verlauft die Bewegung eines Pendelkorpers unter den Vereinfachungen:keine Berucksichtigung des Luftwiderstandes, Masse des Fadens wird vernachlas-sigt, Pendelkorper ist Massepunkt (mathematisches Pendel).

Bekanntlich gilt nach NEWTON F = ma.Fur kleine Auslenkungen α(t) ergibt sich aus dem Newtonschen Grundgesetzfolgende Gleichung fur das mathematische Pendel (Krafteparallelogramm)

m g sinα(t) = −m l α′′(t),

l Fadenlange,l α(t) Weg, den das Pendel zurucklegt = Lange des Kreisbogens,a = l α′′(t) Beschleunigung des Pendelkorpers,g Erdbeschleunigung,mg Gewichtskraft,F und α haben entgegengesetzte Richtungen.

Fur kleine Winkel ist sinα ≈ α, so daß eine einfache lineare DGL 2. Ordnungentsteht.

α′′(t) +g

lα(t) = 0.

U

+?

+

M

l

m g

F

α

α

”– Richtung“

”+ Richtung“

P

s

t

*

(a) Allgemeine Losung der DGL

Hinweis: Finden von Basislosungen mittels Ansatz eλt.

Die beiden unabhangigen Losungen sind sin(√

glt), cos

(√glt).

Die Losungsmenge ist ein Vektorraum gebildet mittels der Linearkombi-nation

α(t) = c1 sin

(√g

lt

)

+ c2 cos

(√g

lt

)

, c1,2 ∈ R.

16

(b) Berucksichtigung von Anfangsbedingungen

Es sollen 2 Varianten unterschieden werden.

• Pendel wird nach Auslenkung um einen Winkel α0 losgelassen.

Wenn zur Zeit t = 0 die AB α(0) = α0 und α′(0) = 0 bekannt sind,ergeben sich die Koeffizienten zu

c1 = 0, c2 = α0.

Daraus folgt

α(t) = α0 cos

(√g

lt

)

.

• Pendel wird aus der Lage α(0) = 0 angestoßen mit der Geschwindig-keit v0 = lα′(0).Wenn zur Zeit t = 0 die AB α(0) = 0 und α′(0) bekannt sind, ergebensich die Koeffizienten zu

c2 = 0, c1 =

√

l

gα′(0).

Daraus folgt

α(t) =

√

l

gα′(0) sin

(√g

lt

)

.

(c) Losung der DGL 2.Ordnung mittels Einschrittverfahren

Das AWP

α′′(t) = −glα(t), α(0), α′(0) gegeben, t ∈ [0, T ]

kann man mittels α′(t) = β(t) auf ein System von zwei DGL 1.Ordnungtransformieren.

α′(t) = β(t), α(0) gegeben

β ′(t) = α′′(t) = −glα(t), β(0) = α′(0) gegeben.

Nunmehr kann man auf das allgemein notierte DGL-System

x′ = f(t, x, y), x(0) gegeben

y′ = g(t, x, y), y(0) gegeben

das numerische Verfahren anwenden.

xn+1 = xn + hf(tn, xn, yn), x0 = x(0), n = 0, 1, ..., tn = nh, h > 0

yn+1 = yn + hg(tn, xn, yn), y0 = y(0).

Vergleiche die Naherungslosung mit der exakten.

17

(d) 2 Varianten der graphischen Darstellung

• Vergleichende Graphik fur Naherung α0, α1, α2, ... (Polygonzug) undexakte Losung α(t).

Illustration am Beispiel α(t) = α0 cos(√

glt).

6

-0

α0

t

• Vergleichende Graphik fur Phasenkurven- der Naherung (α0, α

′0), (α1, α

′1), (α2, α

′2), ..., α′ = β, (Polygonzug)

- und exakten Losung (α(t), α′(t)).Illustration am Beispiel

α(t) = α0 cos

(√g

lt

)

,

α′(t) = −α0

√g

lsin

(√g

lt

)

,

√g

l= 2, α0 = 1.

-

6

α(t)αn

α′n, α

′(t)

-1 0 1

-2

-1

1

2

s

s

s

s

t=0

t=3π/4

t=π/2

t=π/4

18

7. Simulation eines Doppelpendels [Henning, Kutscha]

Wir betrachten ein idealisiertes Doppelpendel.Zwei Punktmassen gleicher Masse m, die an zwei masselosen Stangen der Langel aufgehangt sind, bewegen sich unter dem Einfluß der Schwerkraft mg. DerLuftwiderstand, die Reibung der Lager und andere Storungen werden ver-nachlassigt.Der Winkel α(t) gibt die Lage des inneren Korpers 1 bezuglich der z-Achse an,γ(t) die des außeren Korpers 2 bezuglich der z-Achse an.Ziel ist die Simulation und Beschreibung des Bewegungsablaufes mittels derWinkel α(t) und γ(t) ohne bzw. mit Storgroßen.Die Ausgangsposition sei gegeben durch α(0), γ(0), α(0), γ(0).

6

-s

u

u

$%

$%

α

γ

0 x

m

m

zα(0) = πγ(0) = πα(0) = 0γ(0) = 0

?

g

Ausgangsposition

6

-s

uu

z

0 x

m

ml

l

α

γ

>

6

Idealisiertes Doppelpendel

(a) Charakterisierung der Geometrie, Rotation und Schwingung sowie DGL

• Geometrie, Lage des inneren und außeren Korpers im R3

r1 = l

sinα0

− cosα

, r2 = l

sinα + sin γ0

− cosα− cos γ

• Energieerhaltungssatz

E = α2 +1

2γ2 + αγ cos(γ − α)− 2 cosα− cos γ

19

• Drallsatz

dLA

dt=

02 sinα + sin γ

0

• System von Differentialgleichungen 2.Ordnung

α =−2 sinα+ γ2 sin(γ − α) + cos(γ − α)[sin γ + α2 sin(γ − α)]

2− cos2(γ − α)

γ =−2 sin γ − 2α2 sin(γ − α) + cos(γ − α)[2 sinα− γ2 sin(γ − α)]

2− cos2(γ − α)

(b) Man lose nun die DGL 2.Ordnung mittels Polygonzugverfahren (EULER-CAUCHY-Verfahren), das fur die DGL y = f(t, y, y) mit Anfangsbedin-gungen die Gestalt hat

yn+1 = yn +∆t yn, n = 0, 1, ..., N, tn = n∆t, y0, y0 gegeben,

yn+1 = yn +∆t yn = yn +∆t f(t, yn, yn),

Schrittweite ∆t > 0 und Anzahl der Schritte N (Simulationsdauer).

Fur obiges System arbeitet das PZV gemaß der Reihenfolge

- Berechnung von α, γ anhand der DGL zum Zeitpunkt tn- αn+1 = αn +∆t αn,γn+1 = γn +∆t γn,

- αn+1 = αn +∆t αn,γn+1 = γn +∆t γn.

(c) Tabellarische Ausgabe der Ergebnisse, eventuell auch als File

(d) Darstellung der berechneten Bewegung in verschiedenen Plots

- Bewegungsablauf des Doppelpendels (der beiden Stangen) im(x, z)-Koordinatensystem (Graphik wie oben),

- Darstellung der Winkelfunktionen α(t), γ(t) entsprechendder Simulationsdauer. Dabei sollen Winkelwerte, die außerhalb von [0, 2π]liegen, mittels Modulo-Rechnung auf dieses Intervall transformiertwerden.

Die nachfolgende qualitative Abbildung entstand mit den ABα(0) = π, α(0) = 0, γ(0) = π, γ(0) > 0.

20

-

6

0 10 20 t

Veranderung des Winkels γ(t)

γ

2π

π

1

(e) Man kann in die Betrachtungen noch eine kleine Storung s einbeziehen(Masse auf der x-Achse, die eine geringe Beschleunigung verursacht).Großenordnungen fur s sind 10−5...10−10.In den Beziehungen wirkt sich diese wie folgt aus.

Energieerhaltungssatz

E = α2 +1

2γ2 + αγ cos(γ − α)− 2 cosα− cos γ + s (2 sinα− sin γ)

Drallsatz

dLA

dt=

02 sinα+ sin γ

0

− s

02 cosα+ cos γ

0

Gestortes System von Differentialgleichungen 2.Ordnung


2− cos2(γ − α)

+ scos γ cos(γ − α)− 2 cosα

2− cos2(γ − α)

γ =−2 sin γ − 2α2 sin(γ − α) + cos(γ − α)[2 sinα− γ2 sin(γ − α)]

2− cos2(γ − α)

+ s2 cos(γ − α)(cosα− 1)

2− cos2(γ − α)

Analog kann man die gestorten Winkelfunktionen mittels PZV berechnen,die Ergebnisse tabellarisch ausgeben sowie den Bewegungsablauf graphischdarstellen.

21

2.5 Mathematik

1. Bestimmung von Kurven mittels Abstandsbedingung I

Gesucht sind alle Kurven mit folgenden Eigenschaften:

Die Tangente an einen Kurvenpunkt P schneide die x-Achse in T . Wie lautetdie Gleichung der Kurve, wenn der Abstand zwischen P und T gleich demAbstand zwischen T und dem Koordinatenursprung O ist fur alle Punkte Pder Kurve?

Aus der Abstandsbedingung |OT |2 = |TP |2 und der Tangentengleichungy − y(x0) = y′(x0)(x− x0) am Punkt P (x0, y0) einer Kurve y = y(x) folgt dieDGL

y′ =2xy

x2 − y2.

-

6

u u

u

x

y

O T (xt, yt)

P (x0, y0)

x0

y0

t


Es ist eine homogene DGL, die unter Nutzung der Substitution u =y

xsowie mittels Partialbruchzerlegung gelost werden kann.Mit y′ = u′x+ u erhalt man

y′ =2 yx

1− ( yx)2

u′x+ u =2u

1− u2

u′x =u(1 + u2)

1− u2

1

x=

1− u2

u(1 + u2)u′ =

u′

u− 2uu′

1 + u2

ln |u| − ln |1 + u2| = ln |x|+ c

22

∣∣∣∣

u

1 + u2

∣∣∣∣

= c1|x|, c1 ∈ R+.

Mit u = y/x erfolgt nach Umformung

y = c1(y2 + x2)

(

y − c22

)2

+ x2 =(c22

)2

, c2 = 1/c1.

Das sind offenbar Kreisgleichungen. Die Mittelpunkte der Kreise liegenauf der y-Achse.

6

-s s

sr r

sO T x

P

c22

y

S

α

y = |x|

c2

(b) Man lose nun die DGL y′ = 2xyx2−y2

= f(x, y) mittels Polygonzugverfahren

(EULER-CAUCHY-Verfahren)

yn+1 = yn + hf(xn, yn), n = 0, 1, ..., xn = x0 + nh, y0 = y(x0),

Schrittweite h > 0 gegeben.

Jedoch sind zunachst einige Einschrankungen und Besonderheiten zu beruck-sichtigen.

• Im kartesischen Koordinatensystem kann man mit der DGL wegender Eindeutigkeit der Losungskurve maximal den unteren oder obe-ren Halbkreis beschreiben. Ansonsten ist die DGL zu transformieren,indem die Koordinaten in Abhangigkeit vom Winkel α geschriebenwerden, z.B. x(α) = c2

2sinα, y(α) = c2

2(1− cosα), α ≥ 0.

• Die AB kann nicht (x0, 0) sein.Dann ist y′(x0) = f(x0, 0) = 0 fur x0 6= 0 oder man wurde y′(0) =f(0, 0) per Definition zu Null setzen. Das PZV yn+1 = yn(1+h

2xn

x2n−y2n

)liefert in beiden Fallen die bekannte Nullosung der DGL(Grenzfall = unendlicher Kreis).

23

• Als AB ist y(x0) = |x0|, x0 6= 0, nicht zugelassen, weil dann dieFunktion f(x, y) fur den Startwert nicht definiert ist. Diese AB liegendort, wo die Tangenten an den Kreis vertikal sind.

• Man wahle also geeignete AB S(x0, y0), berechne mittels PZV dieNaherungslosung und kontrolliere ihren Verlauf.Mache einen graphischen Vergleich der exakten Losung (Kreisbogen)mit dem Losungsverlauf des Naherungsverfahrens (Polygonzug).Teste verschiedene Schrittweiten h.

2. Bestimmung von Kurven mittels Abstandsbedingung II

Gesucht sind alle Kurven mit folgender Eigenschaft:

Der Schnittpunkt S der Kurvennormalen n in einem Kurvenpunkt P mit derx-Achse soll von der Abszisse von P stets den Abstand 2 haben.

Aus der Normalengleichung y − y(x0) = − 1

y′(x0)(x− x0) im Punkt P (x0, y0)

einer Kurve y = y(x), ihrem Schnittpunkt S(xs, 0) mit der x-Achse sowie derAbstandsbedingung xs − x0 = 2 folgt die DGL

y y′ = 2.

-

6

u

u

y

0 x0 S(xs, 0) x

Py0

2

n


Es gilt

1

2(y2)′ = 2

y2 = 4x+ c

y = 2√x+ c1, x ≥ −c1.

24

(b) Charakterisiere das Richtungsfeld der DGL y′ = 2y.

In jedem Punkt P (x, y) des Gebietes Ω ∈ R2 kann man die Steigung y′

durch Eintragen eines kleinen”Strichs“ fur die Tangentenrichtung andeu-

ten. Ein solcher Strich heißt Linienelement oderRichtungselement (x, y, y′)der Losungskurve. Die Gesamtheit aller Linienelemente heißtRichtungsfeld.Was passiert, wenn y = 0?

(c) Bestimme die orthogonalen Trajektorien der Kurvenschar.

Definition: Alle Kurven, die in jedem Punkt eine vorgegebene Kurvenscharorthogonal schneiden, heißen orthogonale Trajektorien.

Gegeben sei die explizite DGL 1.Ordnung y′ = f(x, y).

Fur ihre Losungsschar bestimmt man die orthogonale Kurvenschar, indemman lediglich y′ durch − 1

y′ersetzt. Die Vektoren (1, y′) und (1,−1/y′) sind

namlich orthogonal.Somit lautet deren DGL y′

y= −1

2, ihre allgemeine Losung ist y = c1 e

−x

2 .

Zeichne ebenfalls das Richtungsfeld der orthogonalen Trajektorien.

(d) Man lose nun die DGL y′ =2

y= f(x, y) sowie die der orthogonalen

Trajektorien y′ = −y2

mittels Polygonzugverfahren (EULER-CAUCHY-

Verfahren)

yn+1 = yn + hf(xn, yn), n = 0, 1, ..., xn = x0 + nh, y0 = y(x0),


Man wahle geeignete AB (x0, y0), berechne mittels PZV die Naherungslosungund kontrolliere ihren Verlauf.Mache einen graphischen Vergleich der exakten Losung mit dem Losungs-verlauf des Naherungsverfahrens (Polygonzug).Teste verschiedene Schrittweiten h.

3. Phasenkurve

Fur die Losung y = y(x) der DGL y′ = f(x, y) ist manchmal auch eine andereDarstellung zu empfehlen. Das trifft zum Beispiel auf Bewegungsdifferential-gleichungen bei Abhangigkeit von der Zeit t zu. Formal kann man x = t setzenund erhalt das DGL-System

x = 1, y = f(t, y).Die Losung (x(t), y(t)) wird Phasenkurve genannt und ihre Darstellung kannim Phasenraum (x, y) in Abhangigkeit von der Zeit t erfolgen.Mit einer gegebenen Anfangsbedingung (y(t0), x(t0)) = (y0, x0) lost man dasAWP. Allgemein lautet das DGL-System x = f1(t, x, y), y = f2(t, x, y).

Der eigentliche Typ von Phasenkurven entsteht bei der DGL 2. Ordnung

25

y′′ = f(t, y, y′), wenn man diese durch die Definition x = y′ = y in das Systemvon zwei Gleichungen

y = x = f1(t, y, x), x = y = f(t, y, x) = f2(t, y, x)

uberfuhrt und somit die Phasenkurve (y, x) = (y, y) erhalt.

Konstruiere fur folgende DGL die Phasenkurven.

a) y′ = 1 + y, y(0) = 1/2.

b) y′′ + y = 0, y(0) = 0, y′(0) = 1, Pendel ohne Reibung.

c) y′′ + ky′ + y = 0, y(0) = y0, y′(0) = y′0, Pendel mit Reibung k > 0,

z.B. k = 103.

d) y′′ − (2 + et)y′ + y = 0, y(0) = 2, y′(0) = −1.Liefert (1,0) als Grenzwert der Phasenkurve eine stationare Losung?

4. Van der Polsche DGL [Arnold]

Normierte DGL 2.Ordnung

y − ε(1− y2)y + y = 0, y(0) = y0, y(0) = y0, t ≥ 0, ε ≥ 0.

Parameter ε legt die Einschwingdauer fest.

(a) Notation als Phasenkurve

y = ε(1− y2)y − y = f(t, y, y), x ≡ y

System von DGL 1.Ordnung

y = x = f1(t, y, x), y(t0) = y0, t0 = 0

x = f(t, y, x) = f2(t, y, x), x(t0) = y0

Phasenkurve: (y, x) = (y, y).

(b) Van der Pol DGL und Losung fur ε = 0

y + y = 0

y(t) = y0 sin t + y0 cos t

= A sin(ωt+ ϕ), harmonische Schwingung

y(t) =√2 sin(t+

π

4), Losung fur AB y0 = y0 = 1

-

6y

1

0

-1

-π4 π 2π t

Amplitude A =√2

Frequenz ω = 1Periode 2π

ω= 2π

26

(c) Losung der DGL y + ω2y = 0, ω > 0, ε = 0

y(t) = A sin(ωt+ ϕ), Ansatz zur Bestimmung von A, ϕ

y(t) = Aω cos(ωt+ ϕ)

y(t) = −Aω2 sin(ωt+ ϕ)

Berucksichtigung der AB ohne Betrachtung von Sonderfallen

y0 = A sin(ωt0 + ϕ)

y0 = Aω cos(ωt+ ϕ)

tan(ωt+ ϕ) =ωy0y0

ϕ = arctanωy0y0

− ωt0

ω2y20 + y20 = A2ω2 sin2(ωt+ ϕ) + A2ω2 cos2(ωt+ ϕ)

A2 = y20 +y20ω2

A =1

ω

√

(ωy0)2 + y20

Losung

y(t) =1

ω

√

(ωy0)2 + y20 sin

(

ω(t− t0) + arctanωy0y0

)

(d) DGL y + ω2y = 0, ω = 1, ε = 0 als DGL-System

x ≡ y

x = y = −y

Somit y = x, y(t0) = y0x = −y, x(t0) = y0

bzw. x1 = x2, x1(t0) = x10 = y0x2 = −x1, x2(t0) = x20 = y0

Phasenkurve: (y, x) = (y, y) = (x1, x2).

Losungsansatz mit Bestimmung von A aus AB ohne Betrachtung vonSonderfallen

x1(t) = A cos(t + ϕ)

x2(t) = −A sin(t+ ϕ)

27

x1(t0) = A cos(t0 + ϕ) = y0

x2(t0) = −A sin(t0 + ϕ) = y0

φ = − arctany0y0

− t0 = −arccoty0y0

− t0

A =√

y20 + y20

y(t) = x1(t) =√

y20 + y20 cos

(

t− t0 − arccoty0y0

)

=√

y20 + y20 cos

(

t− t0 −(π

2− arctan

y0y0

))

=√

y20 + y20 sin

(

t− t0 + arctany0y0

)

Spezialfalle:

1. y0 = A, y0 = 0 → arctan∞ = π2

y(t) = A sin(

t− t0 +π

2

)

= A cos(t− t0)

2. y0 = 1, y0 = 2 → A =√5, ϕ = − arctan 2− t0 ≈ −1.1− t0

y(t) =√5 sin(t− t0 + arctan 0.5) ≈

√5 sin(t− t0 + 0.46)

y(t0) =√5 sin(0.46) ≈ 1

y(t) =√5 cos(t− t0 + arctan 0.5) ≈

√5 cos(t− t0 + 0.46)

y(t0) =√5 cos(0.46) ≈ 2

-

6

-

6x2

x10 A

x2

2

1

0 1√5 x1

x1,2(t0)

s

s

r

1. x1(t) = A cos(t− t0) 2. x1(t) =√5 cos(t− t0 − 1.1)

x2(t) = −A sin(t− t0) x2(t) = −√5 sin(t− t0 − 1.1)

28

(e) Van der Polsche DGL

y = ε(1− y2)y − y = f(t, y, y), x1 ≡ y,

als System von DGL 1.Ordnung

x1 = x2

x2 = −x1 + ε(1− x21)x2

Damit erhalten wir eine besondere Form der allgemeinen Storungsglei-chung

x1 = x2 + εf1(x1, x2)

x2 = −x1 + εf2(x1, x2)

mit entsprechenden Funktionen f1, f2.

(f) Beispiele zum Schwingungsverhalten, qualitative Verlaufe y(t)

- -

6 6

t t

2 2

-2 -2

0 10 20 0 10 20

1. y − (1− y2)y + y = 0 2. y − 5(1− y2)y + y = 0

y(0) = 0, y(0) = −0.05 y(0) = 2, y(0) = 0

Wachst ε > 0 an, so verzerrt sich die Schwingung immer mehr und diePeriode vergroßert sich.

Bei ε ≫ 1 (z.B ε = 100, 1000) beschreibt die DGL Relaxationsschwin-gungen. Dabei wechseln sich Zeiten sehr langsamer Bewegung mit solchensehr schneller Bewegung ab.Der Parameter ε legt auch die Einschwingdauer fest.

29

2.6 Biologie und Populationsdynamik

1. Zweispezimodell - okologisches Gleichgewicht bei Popu-lationen

Dies wird beschrieben durch das Volterrasche System von zwei DGL 1.Ordnung

x = ax(1 − y), x = x(t),

y = −cy(1− x), y = y(t),

bzw. entkoppelt in der Form

xx = x2 + acx2 − cxx+ cx2x− acx3,

yy = y2 + acy2 + ayy − ay2y − acy3.

Lose die DGL fur a = 2, c = 1 mit den AB x(0) = 1, y(0) = 3.

2. Lotka-Volterra-Modell [Braun]

Das originale Lotka-Volterra-Modell

x = x(1 − y),

y = −y(1− x)

ist ein Rauber-Beute-System mit der Beutepopulation/Futtervorrat x(t) undder Rauberpopulation/Verbraucher y(t) und kann auf eine hohere Dimensionverallgemeinert werden.Die 3D-Verallgemeinerung hat die Form

x1 = x1 − x1x2 + cx21 − ax3x21 = x1(1− x2) + x21(c− ax3),

x2 = −x2 + x1x2 = −x2(1− x1),

x3 = −bx3 + ax1x23 = −x3(b− ax1x3)

mit a, b, c ≥ 0 als Bifurkationsparameter.

3. Konkurrenz zweier Arten, die sich aus derselben beschrank-ten Resource ernahren [Braun]

a) u = u[−a1 + c1(1− b1u− b2v)], u(0) = u0, c1 > a1

v = v[−a2 + c2(1− b1u− b2v)], v(0) = v0, c2 > a2

b) u = u[b1 − λ1(h1u+ h2v)], u(0) = u0

v = v[b2 − λ2(h1u+ h2v)], v(0) = v0

30

4. Jager-Beute-Modell [Kubicek, Marek, Braun]

a) u = u

[

R(

1− u

K

)

− Rv

u+D

]

, u(0) = u0

v = v[

−S(

1− γv

u

)]

, v(0) = v0

b) x = ax− bxy − ex2, x(0) = x0

y = −cy + dxy − fy2, v(0) = y0

c) x = ax− b√x y, x(0) = x0 > 0

y =

d√x y, falls x 6= 0

−cy, falls x = 0, y(0) = y0 > 0

5. Ameisen-Kampf [Braun]

S = 110S − 1

20SN, S(t0) = S0 > 0

N = 1100N − 1

100N2 − 1

100SN, N(t0) = N0 > 0

a) Man zeige, daß die Geraden N = 2 und N + S = 1 den ersten Quadrantenin drei Gebiete aufteilen, in denen dS/dt und dN/dt feste Vorzeichen haben.

-

6

0 1 S

1

2

NIII S < 0, N < 0

II S > 0, N < 0

I

S > 0N > 0

b) Man verifiziere, daß jede Losung S(t), N(t), die in den Gebieten I oder IIIstartet, schließlich nach II geht.

c) Man beweise, daß jede im Gebiet II beginnende Losung S(t), N(t) fur allespateren Zeitpunkte dort bleibt.

d) Man schließe aus c), daß fur alle Losungen S(t), N(t) mit positivenS(t0), N(t0) gilt limt→∞ S(t) = ∞, limt→∞N(t) endlich mit Grenzwert≤ 2.

e) Man beachte, daß ein t0 existiert, so daß dN/dt ≤ −N fur t ≥ t0 gilt. Manschließe aus dieser Ungleichung auf limt→∞N(t) = 0.

31

6. Veranderungen bei Populationen - Wachstumsmodell

Unter der Annahme, daß es eine maximale Population m gibt, die z.B. ausGrunden des Lebensraumes, nicht uberschritten werden kann, sei die Bestands-anderung der Population proportional sowohl zum Bestand als auch zur jewei-ligen Abweichung vom moglichen Hochststand. Dies ergibt die DGL

y′ = f(t, y) = p y(m−y), y = y(t) ≥ 0, AB y(0), Proportionalitatsfaktor p > 0.

Sie wird auch als Wachstumsgleichung oder logistische Gleichung bezeichnet.

(a) Losung des Anfangswertproblems

Hinweis: Trennung der Variablen, Idee der Partialbruchzerlegung.

p =1

y(m− y)y′ =

[1

my+

1

m(m− y)

]

y′

mp =y′

y+

y′

m− y∫

mpdt =

∫y′

ydt+

∫y′

m− ydt

mp t+ c = ln y − ln (m− y) = lny

m− y

c1emp t =

y

m− y, c1 = ec > 0

Somit haben wir eine Bedingung fur das Wachstum, die nach y aufzulosenist.

y(t) =c1me

mpt

1 + c1 empt, c1 ∈ R

+

Unter Berucksichtigung der Anfangsbedingung ergibt sich folgende Glei-chung fur c1.

y(0) =c1m

1 + c1und somit c1 =

y(0)

m− y(0)

Die Wachstumsgleichung hat die endgultige Form

y(t) =y(0)mempt

m− y(0) + y(0) empt.

(b) Eigenschaften der DGL und Losung

y(t) =y(0)mempt

m+ y(0)(empt − 1), t ∈ [0, T ], AB y(0) ∈ (0, m), p > 0

32

-

6

0 t

m

y(0)sy(t)

• Uberprufe folgende Aussagen zur Losung.

y(0) ≤ y(t) < m, m = limt→∞

y(t)

y(t1) < y(t2) fur t1 < t2

y′(t) > 0, y′(0) = p (m− y(0)) y(0)

• Was passiert, wenn die AB y(0) = 0 bzw. y(0) = m ist?

• Charakterisiere das Richtungsfeld der DGL.

In jedem Punkt P (t, y) des Streifens Ω = [0, T ] × [0, m] ∈ R2 kann

man die Steigung y′ durch Eintragen eines kleinen”Strichs“ fur die

Tangentenrichtung andeuten. Ein solcher Strich heißt Linienelementoder Richtungselement (t, y, y′) der Losungskurve. Die Gesamtheitaller Linienelemente heißt Richtungsfeld.Nimm zwecks Vereinfachung die Parameter p = m = 1.Betrachte auch das Richtungsfeld außerhalb des Streifens.

(c) Lose das AWP im Bereich mit zulassigem Startpunkt y(0) mittels Poly-gonzugverfahren (EULER-CAUCHY-Verfahren)

yn+1 = yn + hf(tn, yn), n = 0, 1, ..., tn = nh, y0 = y(0),


Mache einen graphischen Vergleich der exakten Losung mit dem Losungs-verlauf des Naherungsverfahrens (Polygonzug). Teste verschiedene Schritt-weiten h.

Wie ist h zu wahlen, damit z.B. analog zu den Eigenschaften der exaktenLosung die Bedingungen

yn < yn+1 < ... < m

erfullt werden.Rechne das PZV auch mit Anfangsbedingungen, die außerhalb des Strei-fens liegen.

33

7. Modellierung eines Okosystems

Das typische Rauber-Beute-Modell kann auf der Grundlage von gewohnlichenDifferentialgleichungen beschrieben werden. Dazu nimmt man an, daß- die raumliche Verteilung der Spezies keine Rolle spielt,- die Modellierung mit reellen Zahlen aufgrund einer hohen Anzahl vonIndividuen sinnvoll ist,

- eine Normierung der Anzahlen auf 1 vorgenommen wird.

Bezeichnet man die Menge Gras, die Zahl der Hasen und die Zahl der Fuchsemit g, h bzw. f , so kann man folgendes Modell (System von DGL) aufstellen.

g′(t) = 1− h(t) g(t),

h′(t) = h(t)(g(t)− f(t)),

f ′(t) = f(t) h(t)− c1f(t)− c2√

f(t), c1 + c2 = c = const > 0, ci ≥ 0.

Mogliche Anfangsbedingungen (AB) sind g(0) = h(0) = 1, f(0) = 0.1.Die erste Gleichung beschreibt das Wachstum der Grasmenge und das “Gefressen-werden von Gras“ mit einer Rate, die proportional zum Produkt aus Hasen undGras ist. In der zweiten Gleichung fur Hasen ist die Wachstumsrate ebenfallsmit h(t)g(t) gegeben. Zusatzlich geht die Zahl der Fuchse “hasenmindernd“ mit−h(t)f(t) ein. In der Gleichung der Fuchse ist ihr Wachstum von den vorhande-nen Hasen abhangig, wahrend der negative Term das Jagerverhalten simulierensoll. Der Abschuß erfolgt proportional zur Anzahl der Fuchse.

(a) SimulationslaufDas System tendiert zu Schwingungen.Dies werden jedoch mit der Zeitherausgedampft, und es stellt sich ein Gleichgewicht ein.

Beispiel 1: c = c1 = 0.8, c2 = 0

0 5 10 15 20 25 30 0

0.5

1

1.5

2

2.5 Simulation : Oeko-System

t

Gras Fuchs

Hase

34

Beispiel 2: c = 0.8, c1 = 0.3, c2 = 0.5

Wegen 1 >√

f(t) > f(t) am Anfang haben wir eine hohere Abschußquo-te der Fuchse als im ersten Beispiel. Die Zahl der Fuchse wird so starkdezimiert, so daß es zu einer richtigen Hasenplage kommt. Die Folge dazuist wieder ein starker Ruckgang der Grasmenge. Insgesamt sind die Aus-schlage großer.

0 5 10 15 20 25 30 0

0.5

1

1.5

2

2.5

3

3.5

4 Simulation : Oeko-System

t

Gras

Fuchs

Hase

(b) Stationares GleichgewichtBei hinreichend großer Zeit stellt sich fur die Anzahlen der Populationenein Gleichgewicht ein. Berechne die Gleichgewichtsgroßeng∞ = g(t = ∞), h∞ und f∞ aus den Bedingungen des stationaren Ver-haltens g′ = h′ = f ′ = 0.

(c) NumerikDie numerische Losung des Systems der DGL 1.Ordnung

y′ = F (t, y), y : R → R3,

erfolgt mittels Polygonzugverfahren (EULER-CAUCHY-Verfahren)

yn+1 = yn +∆t F (tn, yn), n = 0, 1, ..., N,

tn = n∆t, Schrittweite ∆t > 0 gegeben,

AB y0 = y(0) gegeben.

Man wahle geeignete AB y0, berechne mittels PZV die Naherungslosungenund kontrolliere ihren Verlauf.Stelle die Losungen graphisch dar. Teste verschiedene Schrittweiten ∆t.

Die Abbildungen sind mit Mathematica erstellt.

35

2.7 Medizin

1. Modell der Diabeteserkennung [Braun]

u = F1(u, v) + J(t), u(t) Konzentration der Glukose im Blut

v = F2(u, v), v(t) Insulinkonzentration

J(t) externe Rate des Anstiegs der Blutzuckerkonzentration

2. Epidemiologie [Braun]

Einfaches Modell zur zeitlichen Verbreitung von Epidemien ohne Betrachtungvon Migrationen der Gesamtpopulation.Fragen: Entsteht eine Epidemie oder nicht? Wann stirbt sie aus?

y1 = −ry1y2, y1 aufnahmefahige Population

y2 = ry1y2 − γy2, y2 infizierte Population

y3 = γy2, y3 resistente bzw. gestorbene Population

3. Modell eines Immunsystems [Kubicek, Marek]

Das Modell der Abwehr des Immunsystems gegen Krebs beschreibt die zeitlicheund raumliche Entwicklung der Zellular-Population X(t, r), M0(t, r), M1(t, r).

X(t, r) Konzentration der normalen ZellenM0,1(t, r) Konzentrationen der zytotoxischen Zellen

4. Zellwachstum

Die DGL fur explosionsartiges Wachstum von Zellen habe die Formy′ = y2, y(0) = c > 0, t ≥ 0.

Losungy′

y2= −

(1

y

)′

= 1, y 6= 0

1

y= −t + C, AB beachten

y(t) =1

−t + 1/c

-

6

0 1 t

1

2

y

Bsp.: c = 1, y(t) = 1−t+1

Analytische Eigenschaftender Losung beachten,hier z.B. Polstelle bei t = 1.

36

5. Tumorwachstum [Braun]

Die einfache Dynamik des Zellwachstums kann beschrieben werden durch dieDGL

y′ = λy, y(t0) = y0, λ > 0 Wachstumsrate.

Gemaß ihrer Losung y(t) = y0 exp(λ(t− t0)) wachsen die sich teilenden Zellenexponentiell an, wobei sich das Zellvolumen y(t) wahrend eines jeden Zeitin-tervalls der Lange ln 2/λ verdoppelt.Empirische Untersuchungen haben gezeigt, daß die Zunahme des Zellvolumenseher der Bedingung

y(t) = y0 eλ

α(1− e−αt), λ, α > 0

entspricht (Gompertzsche Beziehung), die fur t→ ∞ den Grenzwert y0 exp(λ/α)hat.Die dazugehorige DGL ist

y′ = λe−αty,

die zweierlei interpretiert werden kann.

a) y′ = (λe−αt) y

Die Hemmung des Zellwachstums hat ihre Ursache in der mit der Zeit fal-lenden Wachstumsrate λ exp(−αt).

b) y′ = λ (e−αty)

Die durchschnittlich zur Zellbildung notige Zeit bleibt konstant, aber dieWachstumshemmung basiert auf dem Verlust sich reproduzierender Zellenexp(−αt)y.

37

3 Implementierung der Modelle

1. Allgemeine Hinweise

Die zahlreichen Modelle eignen sich auch fur Belegaufgaben und/oder ein Com-puterpraktikum. Dabei konnen sowohl Untersuchungen zum Problem als auchdie entsprechende Algorithmen- und Programmentwicklung unter Anwendungder Kenntnisse uber Programmierkonzepte wie Unterprogrammtechnik, Modul-konzept, Rekursion, Programm-Tools, Units, statische und dynamische Daten-strukturen, graphische Bausteine, objektorientierte Programmierung, Entwurf-soberflache Turbo Vision, Menutechnik, Klassenbibliotheken, Einbinden vonProgrammbausteinen, Hilfen u.a. erfolgen.Bei Vergleich von verschiedenen numerischen Methoden, wobei oft mit einfa-chen Verfahren begonnen wird, ist eine Aufwandsananlyse bzgl. Arithmetik,Zeit und Speicherbedarf moglich.Empfehlenswert ist eine Aufwartskompatibilitat der Programme, d.h. zusatzli-che gleichartige Aufgabenstellungen sollten ohne großeren Aufwand integrierbarsein.

Bezuglich des Programmentwurfs orientiere man sich soweit wie moglich an derStruktur des Problems. Dabei sind 3 Entwurfsmethoden hervorzuheben.

• Funktionsorientierter Entwurf, auch schrittweise Verfeinerung oder algo-rithmische Dekomposition genannt.Ausgangspunkt sind Funktionen bzw. Operationen des Systems. Imple-mentierungsdetails, insbesondere Einzelheiten der Realisierung auf demRechner, werden auf die

”untersten“ Module verschoben.

• Datenorientierter Entwurf.Ausgangspunkt sind Operanden, und man orientiert sich an den Daten-strukturen.

• Objektorientierter Entwurf.Wichtiges Element ist dabei die Verknupfung von Operanden und Opera-tionen zu Objekten.

2. Hinweise zur Menutechnik

Die Struktur und der Inhalt der Programmenus konnen mehr oder wenigerelegant gestaltet werden. Sie werden auf dem Textbildschirm oft bescheidenerausfallen als auf graphischen Oberflachen.Der Aufbau einer hierarchischen Menustruktur, manchmal noch mit zahlrei-chen Querverweisen, fuhrt meistens zu einer erheblichen Vergroßerung des Pro-grammkodes. Wichtig ist jedoch, daß der Softwareanwender zielgerichtet auf dieeigentliche Losung des Problems gefuhrt wird und dabei die notwendigen Hil-festellungen erhalten kann. Oft richtet man Menu- und Statuszeilen, Anzeige-,Dialog-, Schalter- Hilfe- und Infofenster ein.

38

Empfohlen wird zum Beispiel, folgende Teilaufgaben des Algorithmus zu einerEinheit zusammenzufassen.

• Initialisierung/Eingabe - Init

- Standardvorgaben ubernehmen- Tastatureingabe- Dateieingabe- Auswahl von Funktionen im Quelltext- Funktionseingabe zur Laufzeit- Schalterwahl, Tastenwahl, Hotkeys- Korrektur und Modifikation der Daten

• Verfahrens- und Methodenauswahl - Proc

- Klassifizierung der Verfahren- Eingabe spezieller Verfahrensparameter

• Rechnen - Calc

- Speichern von Zwischenergebnissen und Ergebnissen- Ablaufsteuerung- Unterbrechungsmoglichkeiten mit Datenrettung- Aktivitatsanzeige

• Anzeigen - Show

- Ergebnisdarstellung- Graphische Komponenten- Ausblenden von Teilen der Ergebnisse- Wiederholung der Anzeige/Ausgabe

• Druck - Print

- Druckerauswahl und -einstellungen

• Dateiarbeit - File

- Festlegung von Dateitypen

• Einstellungen - Options

- Einstellungen fur Text- und Graphikbildschirm- Cursor- und Mausbedienung- Rechnerspezifische Einstellungen

• Wiederholungen - Restart

- Neues Beispiel rechnen- Anzeige der zuletzt abgespeicherten Ergebnisse

• Ruckkehr zu anderen Menus/Programmende - Exit

39

3. Illustration an einem Beispiel

Bezuglich der Programmoberflache ist das nachfolgende Beispiel als Anregungzu sehen. Weitere Moglichkeiten fur Menutechniken kann der Leser dem Skript“Praktikumsaufgaben fur die Fachgebiete Algorithmen und Programmierung,

Grundlagen der Informatik und Informatik und Programmierung“ [Neundorf,

Walther] entnehmen.

Simulation des Doppelpendels aus Abschnitt 2.4.

Autor: C.Reimann, Matrikel 96 Studiengang Technische Physik.Die Implementierung des Pendels ohne/mit Storung erfolgte mittels einem mo-dularen Konzept mit objektorientierten Komponenten in der Programmierspra-che Turbo Pascal. Es lauft auf DOS-Ebene von 80387-PC und hoher. Die hard-waremaßige VESA-Unterstutzung ist empfohlen, ansonsten ist das ProgrammUNIVBE.EXE zu laden.

Teile des Programmsystems sind:

- DOPPEL.PAS Hauptprogramm (Rahmen)- VESA.PAS Routinen zur Programmierung des VESA-Standards

Include-Files

- TYPE.INC Globale Typdeklarationen- CONST.INC Globale Konstanten- VAR.INC Globale Variablen- FORWARD.INC Forward-Deklarationen- INIT.INC Grundinitialisierungen- UTILS.INC Utilities- GRAPH.INC Graphische Oberflachenroutinen- AUSGABE.INC Anzeigeroutinen- PROG.INC Tastenabfrage und Programmkern- DIFF.INC Beschreibung der Differentialgleichungen- CALC.INC DGL-Loser (PZV)- FILE.INC Filehandlingroutinen

Im Programm gibt es mehrere Menus. Die Auswahl der Menupunkte erfolgtmittels Tasten. Nach Eingabe der Problemparameter sowie Berechnung undSpeicherung der Losung kann z.B. mittels Menupunkt F5 Show folgende gra-phische Losungsdarstellung aufgerufen werden.

40

Graph α γ α′ γ′

-

6

2π

π

0 T ime

Picturebox

Funktionen Startwerte

Aktuelle Zeit: 40 secSpeed: 500

Darstellung der Winkelfunktionen α(t), γ(t), α′(t), γ′(t), 0 ≤ t ≤ 40

Tansformation der Winkelwerte auf das Intervall [0, 2π]

c css

ss

F1 α AN/AUS N.α = 180.0

F2 α′ AN/AUS N.α′ = 0.000

F3 γ AN/AUS N.γ = 180.0

F4 γ′ AN/AUS N.γ′ = 0.001

F5 Weg/Stoer. S.α = 180.0

F6 schneller S.α′ = 0.000

F7 langsamer S.γ = 180.0

F8 Restart S.γ′ = 0.001

F9 Exit ∆t = 2000µs

Das System von DGLn 2.Ordnung fur die zeitabhangigen Winkelfunktionen ist


2− cos2(γ − α)

+ scos γ cos(γ − α)− 2 cosα

2− cos2(γ − α)

γ =−2 sin γ − 2α2 sin(γ − α) + cos(γ − α)[2 sin α− γ2 sin(γ − α)]

2− cos2(γ − α)

+ s2 cos(γ − α)(cosα− 1)

2− cos2(γ − α), s Storungsparameter

Die numerische Losung des Anfangswertproblems wird mittels Polygonzug-verfahren (EULER-CAUCHY-Verfahren) sowohl fur den ungestorten als auchgestorten Fall (N.α(t),... bzw S.α(t),...) durchgefuhrt.

Die Darstellung der berechneten Bewegung erfolgt in verschiedenen Plots

- als Bewegungsablauf der beiden Stangen im (x, z)-Koordinatensystem,- wahlweise Darstellung der Winkelfunktionen α(t), α′(t), γ(t), γ′(t) .

Die Geschwindigkeit des Ablaufs kann durch Steuerparameter geandert werden.

41

4 Allgemeine Softwarehinweise

Aus der umfangreichen Menge von Software zu DGL sollen hier nur einige Program-me, Programmsysteme sowie strukturierte Modul- und Textbook-Bibliotheken aus-gewahlt werden. Die Auswahl erfolgte unter dem Gesichtspunkt der lokalen Verfugbar-keit.Die vielfaltigen Moglichkeiten in Computeralgebrasystemen seien nur erwahnt.

4.1 Programme

1 S GEBIET

• Kurzcharakteristik zum Programm S GEBIET.PASStabilitatsgebiet von MSV fur gewohnliche DGL 1.Ordnung

– (C) W.Neundorf IfMath TUIlmenau 1995

– MS DOS

– Programmiersprache: Turbo Pascal V.6.0,7.0, Borland Pascal V.7.0

– Aufruf: S GEBIET.EXE

– Dokumentation: Quelltext, selbstdokumentierend

– Freeware

– InhaltGraphik des Stabilitatsgebietseines r-stufigen BDF-Verfahrens oder Mehrschrittverfahrens n-ter Ord-nung fur eine gewohnliche Differentialgleichung 1.Ordnung y′ = f(x, y)

Das MSV hat die allgemeine Form

αryr + αr−1yr−1 + ...+ α0y0 = h(βrfr + βr−1fr−1 + ... + β0f0)

fk = f(xk, yk), y0, y1, ..., yr−1 gegeben

Spezielle einbezogene explizite und implizite MSV sind:

BDF, GEAR

αryr + αr−1yr−1 + ...+ α0y0 = hβrfr, βr = 1, falls normiert

ADAMS-BASHFORTH

αryr = −αr−1yr−1 + h(βr−1fr−1 + ...+ β0f0), βr = 0

ADAMS-MOULTON

αryr = −αr−1yr−1 + h(βrfr + βr−1fr−1 + ...+ β0f0)

42

NYSTRoM

αryr = −αr−2yr−2 + h(βr−1fr−1 + ...+ β0f0), βr = 0

MILNE-SIMPSON

αryr = −αr−2yr−2 + h(βrfr + βr−1fr−1 + ...+ β0f0)

• Kurzcharakteristik der Oberflache

– Menuorientiert, wahlweise Ablaufsteuerung

– Auswahl zwischen vorgegebenen MSV und BDF (meist bis zur Ordnung6), berechneten BDF (maximal bis zur Ordnung 20) und anderen MSVmit Koeffizienteneingabe

– Graphik des Stabilitatsgebiets des MSV mit Auswahl des Gebiets und derSkalierung

– Dateiarbeit, Stabilitatsgebiet (Rand) als Datei verwaltet

– Vergleichende Graphik aller Gebiete einer Klasse von Verfahren, soweitdazu die Dateien vorher erstellt wurden

2 RICHT DG

• Kurzcharakteristik zum Programm RICHT DG.PASRichtungsfeld und AWP fur gewohnliche DGL 1.Ordnung


– MS DOS


– Aufruf: RICHT DG.EXE

– Modul: Unit FOR COM8.PAS fur Funktionsparser f(x, y)


– Freeware

– InhaltGraphik des Richtungsfeldeseiner gewohnlichen Differentialgleichung 1.Ordnung y′ = f(x, y) im auswahl-baren Darstellungsbereich [xa, xa + (j − 1)∆x]× [ya, ya + (n− 1)∆y].Nach Wahl Losung eines AWP im Bereich mit zulassigem Startpunkty(x0) = y0.Druck von ausgewahlten Bereichen des Richtungsfeldes.

43


– Menuorientiert, wahlweise Ablaufsteuerung, benutzerfreundliche Oberflache,Fenstertechnik (4 kleine oder 1 großes Fenster),Ausdunnen des Punktgitters fur Darstellung (fein, normal, grob)

– Demonstrationsvariante

– Auswahl der DGL (Anzahl=20) oder Funktionseingabe und -parser furrechte Seite

– Dateiarbeit, Datei fur (x, y, y′), vorhandene Beispieldateien

– Tabellarische und graphische Ausgaben

– Druck von ausgewahlten Bereichen des Richtungsfeldes

– Losung des AWP mittels klasssischen RKV mit dualer Schrittweitensteue-rung, Auswahl der Parameter, Informationen zum Losungsablauf und Er-gebnissen, Graph der Losung

3 STAB ESV

• Kurzcharakteristik zum Programm STAB ESV.PASStabilitat von ESV bei AWP fur gewohnliche DGL 1.Ordnung


– MS DOS


– Aufruf: STAB ESV.EXE

– Modul: Funktion f(x, y) im Quelltext definiert


– Freeware

– InhaltGraphik des Losungsverhaltens von einfachen ESV furdie gewohnliche Differentialgleichung 1.Ordnung y′ = f(x, y) = y(1− y)(Wachstumsgleichung) im auswahlbaren Argumentbereich [0, xmax].

Test von expliziten Einschrittverfahren- EULER-CAUCHY-Verfahren (PZV)- Klassisches RK-Verfahren (RKV)

yn+1 = yn + hΦ(xn, yn; h), n = 0, 1, 2, ...,

mit y0 ≥ 0 und verschiedenen Schrittweiten h > 0 .

Nach Wahl- Losung eines AWP im Bereich mit zulassigem Startpunkt y(x0) = y0,

44

Losungsverlauf und Stabilitat- Bifurkationsdiagramm und stationare Punkte fur ESV imSchrittweitenbereich 0 < hmin ≤ h ≤ hmax,fur jedes h Berechnung hinreichend vieler Werte yn bis Erreichen desstationaren Zustands, letzte Werte als Datei speicherbar


– Menuorientiert, wahlweise Ablaufsteuerung

– Auswahl zwischen vorgegebenen ESV und Darstellungsvarianten,Auswahl der entsprechenden Parameter

– Graphik der Naherungslosung bzw. des Bifurkationsdiagramms

– Infos zum Losungsablauf und Ergebnissen, Dateiarbeit fur Bifurkations-diagramm

4 DifEqu

• Kurzcharakteristik zum Programm DifEquzur Losung von gewohnlichen DGL

– (C) Geza Makay Ungarn

– Download from the WWW serverhttp : //www.mcs.dundee.ac.uk : 8080/software/index.html#DEssiehe auch http : //www.math.u− szeged.hu/%7Emakay

– MS DOS oder Windows

– Aufruf: difequ.exe

– Dokumentation: Hilfen, selbstdokumentierend

– Freeware

• Inhalt

DifEqu is a programm designed for numerically solving ordinary, functiona-la and partial differential equations, difference equations and do many morethings. It can be used for solving problems arising in mathematics, physics,chemistry, biology. The program ist most usefull for teaching, doing researchand creating simulation. To this end DifEqu has

– A built in programming language in which the user can write a programto be executed by DifEqu.The programming language contains:

∗ Many different variable types (long, float, text, program, file, etc.);

∗ Usual operators (+,−, ∗, / <,>, etc.);

45

∗ User defined variables, vector variables, functions;

∗ Built in functions:- Mathematical functions (Sin, Cos, Log, Exp, Min, Max, etc.);Solving more complicated equations faster then Phaser by Kocak;

- Text handling functions (StrLen, StrMid, Asc, Chr, etc.);- File handling functions (OpenFile, CloseFile, Fseek, FPrint,ReadChar, ReadText, ReadLong, etc.)

- Printing functions (Print, PrintF, AnyToS, Pos, etc.);- Drawing functions (MoveTo, LineTo, PenSize, FgrColor, BgrColor,ChangeFont, etc.);

- Program control functions (Execute, If, While, GetConfig, SetConfig,Abort, etc.);

- Miscellaneous functions (Push, Pop, Clear, SaveValue, etc.)

∗ Descriptive messages and exact location of the error in a programwritten in then programming language;

∗ Program status button to check the values of variables or run a pro-gram by hand.

– Full context sensitive help;

– Many ways of saving your work (position and contents of the windows onthe screen, graphics as Windows BMP file, text as text file);

– Capability of print previewing and printing your work,

– Some examples shipped together with the program.

46

5 DG, DG87

• Kurzcharakteristik zu DG Version 1.1 1994Numerische Losung von GDGL und DGL-Systemen 1.Ordnung

– (C) Frank Schilder IfMath TUI 1994

– POWER-Basic V. 3.0 MS DOS

– Aufruf: DG.EXE StandardversionDG87.EXE mit mathematischen KoprozessorDGL.EXE selbstextrahierende Datei

– Dokumentation: DG.DOC

– Shareware

• Integrierte Verfahren fur eine GDGL 1.Ordnung

– Dormand-Prince-Verfahren der Ordnung 5(4) und 8(7)

– Rosenbrock-Wanner-Verfahren der Ordnung 4(3) und 6(7)

– Gragg-Burlisch-Stoer-Verfahren der maximalen Ordnung 16

– Seulex-Extrapolationsverfahren mittels linear impliziten Euler-

Verfahren der maximalen Ordnung 8


– Integration kleiner und mittlerer AWP

– Menuorientiert, wahlweise Ablaufsteuerung, benutzerfreundliche Oberflache

– On-Line-Hilfesystem

– Funktionseingabe und -parser, Dateiarbeit, vorhandene Beispieldateien

– Tabellarische und graphische Ausgaben

6 APPLI unter Windows

• Kurzcharakteristik zu APPLIGewohnliche Differentialgleichungen

– (C) Falk Nedwal IfMath TUI 1993, 94

– TP fr Windows 1.5 oder BP

– Aufruf: APPLI.EXE

– Dokumentation: README8.TXT

– Shareware

47

• Integrierte Verfahren fur GDGL

– Runge-Kutta-Verfahren der Ordnung 8 (RKV(8))

– Eingebettetes RKV(8,7)

– Rosenbrock-Wanner-Verfahren der Ordnung 5(6) (ROW(5,6))

– Eingebettetes ROW(3,4)

– Seulex-Extrapolationsverfahren mittels linear impliziten Euler-

Verfahren

• Kurzbeschreibung des AWP

– DGL oder DGL-System 1.Ordnung, n ≥ 1, t ∈ [t0, t1], xi = xi(t)

x′1 = f1(t, x1, x2, ..., xn), x1(t0) = x1,0

x′2 = f2(t, x1, x2, ..., xn), x2(t0) = x2,0

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

x′n = fn(t, x1, x2, ..., xn), xn(t0) = xn,0

Variablen konnen auch anders benannt sein, z.B. y′1(x), ...

– Fortsetzung der Integration ab t1 moglich

– Zusatzlich Anzeige von gewunschten ’Observationen’,

z.B. O1(t) = x21 + x22 bei Phasenkurven, O2(t) =n∑

i=1

xi

– Tabellarische Ausgabe

7 DELTA 60

• Kurzcharakteristik zu DELTA 60Dynamische Analyse von gewohnlichen Differentialgleichungen2. und 3.Ordnung, Software zur angegebenen Literatur

– (C) W. Buntig Fak. EI TUI 1992

– TP 6.0 fur MS DOS

– Aufruf: DELTA.EXE

– Dokumentation: DELTA1.DOC

– Literatur: E.Philippow, W.Buntig: Analyse nichtlinearer dynamischerSysteme der Elektrotechnik. Carl Hanser Verlag Munchen 1992.

– Shareware

48

• Hauptmenu

– Wahl der DGL

– Integration der DGL

– Vorbereitung der Fixpunktsuche

– Suche der Fixpunkte

– Fixpunktverfolgung bei Parameteranderung

– Harmonische Analyse

– α− und ω− Invarianzkurven bestimmen

– Leistungsdichtespektrum bestimmen


– Text- und Graphikbildschirme

– Graphik, Speichern und Laden von Bildern

– Wahlweise Ablaufsteuerung

– Auswahlkatalog von vorhandenen DGL, aktuell 16 DGL

z.B. Problem Nr.1 : Van der Polsche DGL

y − ε(1 − y2)y + y = 0, y(0) = y0, y(0) = y0, Parameter ε legt dieEinschwingdauer fest.

Als DGL-System

dx0/dt = x1,dx1/dt− p1(1− x20)x1 + x0 = 0.

8 BVPS unter Unix

• Kurzcharakteristik zu BVPSNumerische Losung von GDGL und DGL-Systemen 1.Ordnung.Objektorientiertes C++ Programm zur numerischen Losung von AWPgewohnlicher DGL-Systeme

x = f(x, t), x(t0) = x0, f : Rn × R → Rn

– Lit.: F.Schilder: Ein portables C++ Programm zu numerischen Inte-gration von Anfangswertproblemen. Preprint M 1/97TUI Marz 1997.

– (C) Frank Schilder IfMath TUI 1996e-mail: [email protected]

– C++ in Unix, auch fur Gnu-C++, Borland-C++ Compiler

49

– Aufruf: bvps dateiname

– Dokumentation: bvps.txt

– Freeware

– siehe auchhttp : //imath.mathematik.tu − ilmenau.de/ ∼ fschild/

• Integrierte Verfahren

– 14 eingebettete explizite ESV hoher Konvergenzordnung• 4 Verfahren von Dormand-Prince

• 4 Verfahren von Fehlberg, (RKF-Verfahren)• 4 Verfahren von Verner

• Verfahren von Merson

• Verfahren von England

– 2 robuste eingebettete linear impliziteRosenbrock-Wanner-Verfahren

insbesondere fur steife DGL

– 3 Extrapolationsverfahren mit Ordnungs- und Schrittweitensteuerung• Limex mittels linear impliziter Mittelpunktsregel• Seulex mittels linear impliziten Euler-Verfahren

• Gragg-Burlisch-Stoer (GBS)


– Integration mittlerer und großerer Systeme, Grenze fur n durch Hardware

– Einfache Bedienbarkeit

– Unkomplizierte plattformunabhangige Oberflache

– Ein- und Ausgaben in ASCII-Format, Dateiarbeit, vorhandene Beispiel-dateien

– Integrierter komfortabler Funktionsparser fur f(x, t)

– Tabellarische und graphische Ausgaben,Ausgabedateien konnen mit Graphikprogramm GnuP lot verarbeitet wer-den.

50

4.2 Programmsysteme

1 MAYA

• Kurzcharakteristik zu MAYA Version 1.0Numerische Methoden

– Lit.: K.Kose, R.Schroder, K.Wieliczek:Numerik sehen und verstehen. Vieweg Braunschweig 1992

– (C) Frank Schilder IfMath TUI 1994

– MS DOS

– Aufruf: MAYA.EXE

– Dokumentation: Selbstdokumentierend

– Shareware

– Inhalt- Funktionen, Interpolation, Konstruktion von Bezier-Polynomen- Ausgleichsrechnung- Chaos bei Differenzengleichungen- Anfangswertaufgaben- Nullstellenaufgaben, Nichtlineare Gleichungssysteme

• Anfangswertaufgaben

– Losungsschar einer Differentialgleichung

– Funktionsweise verschiedener Verfahren

– Stabilitat von Einschrittverfahren

– Vergleich der Verfahren

– Abhangigkeit der Losung von den Anfangswerten

– Zweidimensionale Anfangswertaufgaben (Systeme 1.Ordnung mit 2 DGL)

– Einfluß der Anfangswerte bei zweidimensionalen DGL

• Integrierte Verfahren fur GDGL

– Explizites Euler-Cauchy-Verfahren (Polygonzugverfahren)

– Implizites Euler-Verfahren

– Verbessertes Euler-Verfahren

– Heun-Verfahren

– Trapezregel

– Klassisches Runge-Kutta-Verfahren

51

• Kurzbeschreibung des AWP

– DGL oder DGL-System 1.Ordnung, n = 2, x ∈ [x0, x1]

y′ = f(x, y), y(x0) = y0 bzw.

u′ = f(x, u, v), u(x0) = u0v′ = g(x, u, v), v(x0) = v0

– Verschiedene Graphikkomponenten

2 NUMI2

• Kurzcharakteristik zu NUMI2 Version 1.0 1993Numerische Methoden

– (C) Frank Schilder IfMath TUI, M.Jungierek, S.Thiede, M.Kramer Mer-seburg 1993

– POWER-Basic V. 3.0 MS DOS

– Aufruf: NUMI2.EXE

– Dokumentation: HANDBUCH.TXT

– Shareware

– Inhalt- Gaußscher Algorithmus, LU-Zerlegung, Cholesky-Zerlegung- Iterationsverfahren fur lineare GS, QR-Zerlegung- Hessenberg-Verfahren und Nullstellen des charakteristischen Polynoms- Newton-Verfahren im R und R

n

- Interpolationsverfahren, Spline-Interpolation- Fourier-Interpolation/Entwicklung- Tschebyscheff-Approximation- Differenzenquotient-Verfahren fur Ableitungsbestimmung- Newton-Cotes-Formeln, Gauß-Quadratur- Verbessertes Polygonzugverfahren fur DGL 1.Ordnung- Klassisches Runge-Kutta-Verfahren fur DGL 1.Ordnung- Jacobi-Verfahren fur EWP

• Integrierte Verfahren fur eine GDGL 1.Ordnung

– Verbessertes Polygonzugverfahren

– Klassisches Runge-Kutta-Verfahren


– Menuorientiert, wahlweise Ablaufsteuerung, NUMI-Arbeitsblatt

– Funktionseingabe und -parser, Dateiarbeit

– Graphische Darstellungen

52

4.3 Strukturierte Modulbibliotheken

A. Turbo Pascal Numerical Methods ToolboxBorland International, Inc. 1986, 87

D:\USER\NEUNDORF\MAT123 (lokal)Q:\STUD M93\MAT123 (Novell-Netz MATH1 des IfMath)H:\LANG\TP70\NEUNDORF\MAT123 (RZ-Novell-Netz BLANC)

1. Lineare Algebra

2. Nullstellen von Funktionen

3. Interpolation von Funktionen

4. Approximation von Funktionen

5. Numerische Differentiation und Integration

6. Minimierung und Maximierung von Funktionen

7. Gewohnliche Differentialgleichungen und -systeme

• Anfangswertprobleme (Unit: INITVAL.PAS)

– DGL 1.Ordnung mit Runge-Kutta-VerfahrenUP: InitialCond1stOrder, Programm: RUNGE 1.PAS

– DGL 1.Ordnung mit Runge-Kutta-Fehlberg-VerfahrenUP: RungeKuttaFehlberg, Programm: RKF 1.PAS

– DGL 1.Ordnung mit Adams-Bashforth-Moulton(PECE)-VerfahrenUP: Adams, Programm: ADAMS 1.PAS

– DGL 2.Ordnung mit Runge-KuttaUP: InitialCond2ndOrder, Programm: RUNGE 2.PAS

– DGL n-ter Ordnung mit Runge-KuttaUP: InitialCondition, Programm: RUNGE N.PAS

– DGL-System 1.Ordnung mit Runge-KuttaUP: InitialConditionSystem, Programm: RUNGE S1.PAS

– DGL-System 2.Ordnung mit Runge-KuttaUP: InitialCondition System2, Programm: RUNGE S1.PAS

• Randwertprobleme (Unit: INITVAL.PAS)

– Lineares Schießverfahren fur DGL 2.OrdnungUP: LinearShooting, Programm: LINSHOT2.PAS

– Nichtlineares Schießverfahren fur DGL 2.OrdnungUP: Shooting, Programm: SHOOT2.PAS

• Eigenwertprobleme

8. Partielle Differentialgleichungen und -systeme

9. Integralgleichungen

10. Darstellung spezieller Funktionen

53

B. MATHPAK 87

Source Code for Numerical SoftwareTurbo Pascal 5.0/5.5/6.0 1986-90Precision Plus Software939 Griffith Street London, Ontario Canada N6K 352

Literatur: Mathpak 87 Version 3.0High Performance Software fur den 80x87Benutzer- und ReferenzhandbuchSOS Software Service GmbHAlter Postweg 10, Augsburg

O:\MATPAK (Novell-Netz MATH1 des IfMath)H:\LANG\TP70\NEUNDORF\MATPAK (RZ-Novell-Netz BLANC)

C. TPNUM 3.21 (1993)

Tests, Unterprogramme und Units zur Numerischen MathematikTurbo Pascal, Borland Pascal 7.0Literatur: Dietel,J.(Rechenzentrum der RWTH Aachen):

Formelsammlung zur Numerischen Mathematikmit Turbo Pascal-Programmen (3.Auflage),BI-Wissenschaftsverlag 1991

Q:\STUD M93\FORMELN\TPNUM (Novell-Netz MATH1 des IfMath)O:\FORMELN\TPNUM (Novell-Netz MATH1 des IfMath)H:\LANG\TP70\NEUNDORF\FORMELN\TPNUM (RZ-Novell-Netz BLANC)

∗ Gewohnliche Differentialgleichungen und -systeme

– DGLESV : Einschrittverfahren von Euler-Cauchy, von Heun und vonRunge-Kutta fur eine einzelne Differentialgleichung

– EINBRUKU : Runge-Kutta-Einbettungsformeln mit Schrittweitensteue-rung

– IMPLRUKU : Implizite Runge-Kutta-Verfahren vom Gaußtyp

– Adaptive AnfangswertproblemloserAWP : Einbettungsformel (Runge-Kutta 2. und 3.Ordnung oder Formel

von England 4. und 5.Ordnung)RK FEHL : Runge-Kutta-Fehlberg-Verfahren

– PRAE KO : Pradiktor-Korrektor-Verfahren nach Adams-Bashforth-Moulton

– BUL STOER : Extrapolationsverfahren von Bulirsch-Stoer-Gragg

– RWP : Zuruckfuhrung des Randwertproblems auf ein Anfangswertproblem

54

4.4 Textbook-Bibliotheken

Kockler, N.:

Numerische Algorithmen in Softwaressystemenunter besonderer Bercksichtigung der NAG-Bibliothek.Teubner Stuttgart 1990.

1 Lineare Gleichungssysteme

2 Lineare Optimierung

3 Interpolation und Approximation

4 Nichtlineare Gleichungen

5 Eigenwertprobleme

6 Numerische Integration

7 Anfangswertaufgaben

• Fehler, Konsistenz, Konvergenz, Stabilitat, steife DGL-Systeme

• Einschrittverfahren, Runge-Kutta-VerfahrenNAG-Routine: D02PAF, D02BDF

• Mehrschrittverfahren, Adams-Verfahren, PECENAG-Routine: D02QAF

• Differentiationsverfahren, BDFNAG-Routine: D02QDF

• Anfangswertaufgaben bei IMSL

8 Rand- und Eigenwertprobleme

• DifferenzenverfahrenNAG-Routine: D02GAF

• Schießverfahren, Einfach-Schießen, MehrzielmethodeNAG-Routine: D02HAF

• Randwertaufgaben bei IMSL

9 Partielle Differentialgleichungen

A Die FORTRAN-ProgrammeAnfangswertaufgaben Routine: AWARand- und Eigenwertaufgaben Routine: RWA

B Die NAG-GS Graphik Bibliothek

C PAN - der schnelle Zugriff zu NAG

D Der Aufbau der NAG-Bibliothek

55

Rice, J.R.:

Numerical Methods, Software and Analysis. 2nd Edition.Academic Press Inc. Boston 1993.

1 Mathematics and computer science backround

2 Numerical software

3 Errors, round-off, and stability

4 Models and formulas for numerical computation

5 Interpolation

6 Matrices and linear equations

7 Differentiation and integration

8 Nonlinear equations

9 Ordinary differential equations

• Introduction, difference equations, stability

• Basic method for IVP, EULER, RUNGE-KUTTA, Multistep, PECE

• Polyalgorithms, step size, order and error control

• Systems of ODE

• Software and programms for stiff problemsACM Algorithms, Software Packages IMSL(IVP, Two-Point BVP), DEPAC, EPISODE

• Analysis of methods for ODE, one-step and multistep methods

10 Partial differential equations

11 Approximation of functions and data

12 Software practice, costs, and engineering

13 Software performance evaluation

14 The validation of numerical computations

15 FORTRAN

4.5 Software im Unix-Netz

Im Netz des Instituts fur Mathematik der TU Ilmenau ist durch DM F.Schilder Soft-ware installiert, die unter anderem fur numerische Berechnungen eingesetzt werdenkann. Informationen dazu bietet die WWW-Seite

http : //www.mathematik.tu− ilmenau.de/ ∼ developDort ist auch eine Klassenbibliothek fur Numerik gewohnlicher Differentialgleichun-gen in C++ : fs− classlib sowie die Numerikprogramme- ivps Numerische Approximation zu GDGL,- torus Approximation invarianter 2-Tori.

56

Literatur

[1] Knobloch, J.; Steigenberger, J.: Gewohnliche Differentialgleichungen. Skriptum 100 S.TU Ilmenau 1996.

[2] Schwetlick, H.; Kretzschmar, H.: Numerische Verfahren fur Naturwissenschaftler undIngenieure. Fachbuchverlag Leipzig Koln 1991.

[3] Henning, K.; Kutscha, S.: Informatik im Maschinenbau, 4.Auflage. Springer-VerlagBerlin 1994.

[4] Braun, M.: Differentialgleichungen und ihre Anwendungen. Springer-Verlag Berlin1979, 1991, 1994.

[5] Holmes, P.; Guckenheimer, J.: Nonlinear oscillations, dynamical systems, and bifurca-tion of vector fields. Springer-Verlag 1983, 1990.

[6] Troger, H.; Steindl, A.: Nonlinear stability and bifurcation theory. Springer-Verlag1991.

[7] Seydel, R.: Practical bifurcation and stability analysis. Springer-Verlag New York 1994.

[8] Seydel, R.: From equilibrium to chaos. Elsevier New York 1988.

[9] Goering, H.; Felgenhauer, A.; Lube, G.; Roos, H.-G.; Tobiska, L.: Singularly perturbeddifferential equations. Mathematical Research 13. Akademie-Verlag Berlin 1983.

[10] Lorenz, E.N.: The essence of chaos. Paperback University of Washington Press 1995.

[11] Steeb, W.-H.: Chaos and Fractals. Algorithms and Computation. B.I. Wissenschafts-verlag Mannheim 1992.

[12] Philippow, E.: Grundlagen der Elektrotechnik. Verlag Technik 1988.

[13] Philippow, E.; Buntig, W.G.: Analyse nichtlinearer dynamischer Systeme der ET:Einfuhrung in die numerische Untersuchung einfacher Systeme. Hanser Munchen 1992.

[14] Uhlmann, H.: Grundlagen der elektrischen Modellierung und Simulationstechnik.Akad. Verlagsgesellschaft Geest & Portig. Wiss. Monographien der ET, Bd.5 1977.

[15] Kubicek, M.; Marek, M.: Computational Methods in Bifurcation Theory and Dissipa-tive Structures. Springer Series in Comput. Physics 1983.

[16] Arnold, V.I.: Obyknovennye differentialnye uravnenia. Izd. Nauka Moskva 1971. Dt.Ausgabe : Gewohnliche Differentialgleichungen. Springer-Verlag Berlin 1980.

[17] Hairer, E.; Norsett, S.P.; Wanner, G.: Solving Ordinary Differential Equations, Vol.1,2. Springer-Verlag Berlin 1987, 1991. (Lehrbuchsammlung)

[18] Golub, G.H.; Ortega, J.M.: Wissenschaftliches Rechnen & Differentialgleichungen. Hel-dermann Verlag Lemgo 1995. Hrsg. der dt. Ubersetzung: R.D.Grigorieff.

57

[19] Albrecht, P.: Die numerische Behandlung gewohnlicher Differentialgleichungen.Akademie-Verlag Berlin 1979.

[20] Shampine, L.F.; Gordon, M.K.: Computerlosung gewohnlicher Differentialgleichungen.Vieweg & Sohn Braunschweig 1984.

[21] Fried, I.: Numerical Solution of Differential Equations. Computer Science an AppliedMathematics. Academic Press New York 1979.

[22] Miranker, W.L.: Numerical Methods for Stiff Equations and Singular PertubationProblems. Mathematics and Its Applications 5. D. Reidel Publishing Comp. Dordrecht1981.

[23] Collatz, L.: The Numerical Treatment of Differential Equations, 3rd Edition. Springer-Verlag Berlin 1960.

[24] Abraham, R.H.; Shaw, C.D.: Dynamics - The Geometry of Behavior, Vol. 1, 2, 3(Periodic, Chaotic, Global Behavior) The Visual Mathematics Library. Aerial Press,Inc. Santa Cruz 1983

[25] Parker, T.S.; Chua, L.O.: Practical Numerical Algorithms for Chaotic Systems.Springer-Verlag New York Berlin 1989.

[26] Stuart, A.; Humphries, A.: Dynamical Systems and Numerical Analysis. CambridgeMonographs on Applied and Computational Mathematics 2. Cambridge Univ. Press1996.

[27] Grigorieff, R.D.: Numerik gewohnlicher Differentialgleichungen.Bd. I : Einschrittverfahren. Bd.II : Mehrschrittverfahren. Teubner Stuttgart 1972,1977.

[28] Demailly, J.-P.: Gewohnliche Differentialgleichungen. Theoretische und numerischeAspekte. Vieweg Lehrbuch Angewandte Mathematik Braunschweig 1994.

[29] Strehmel, K.; Weiner, R.: Numerik gewohnlicher Differentialgleichungen. Teubner Stu-dienbucher Mathematik. B.G. Teubner Stuttgart 1995.

[30] Burrage, K.: Parallel and Sequential Methods for Ordinary Differential Equations.Oxford University Press Oxford 1995.

[31] Ainsworth, M.; Levesly, J.; Light, W.A.; Marletta, M.: Theory of Numerics of Ordinaryand Partial Differential Equations. Oxford University Press Oxford 1995.

[32] Korsch, H.J.: Chaos: a program collection for the PC. Springer-Verlag Berlin 1998.

[33] Korsch, H.J.: Chaos: A program collection for the PC. With 250 fig., many numericalexperiments, and two 3 1/2MS-DOS-diskettes. Springer-Verlag Berlin 1994.

[34] Brauning, G.: Gewohnliche Differentialgleichungen. Fachbuchverlag Leipzig 1964.

[35] Werner, H.; Arndt, H.: Gewohnliche Differentialgleichungen. Eine Einfuhrung in Theo-rie und Praxis. Springer-Verlag Berlin 1986.

58

[36] Wenzel, H.: MINOL 7/1,2. Gewohnliche Differentialgleichungen, Band 1, 2. B.G. Teub-ner Verlagsgesellschaft 1973.

[37] Ross, S.L.: Differential equations, 3 Ed. John Wiley & Sons, Inc. New York 1984.

[38] Furlan, P.: Das gelbe Rechenbuch 3 fur Ingenieure, Naturwissenschaftler und Mathe-matiker. Verlag Martina Furlan Dortmund 1997.

[39] Kockler, N.: Numerische Algorithmen in Softwaressystemen : unter besonderer Beruck-sichtigung der NAG-Bibliothek. B.G. Teubner Stuttgart 1990.

[40] Rice, J.R.: Numerical Methods, Software and Analysis. 2nd Edition. Academic PressInc. Boston 1993.

[41] Engeln-Mullges, G.; Reutter, F.: Formelsammlung zur Numerischen Mathematik mitFORTRAN 77- bzw. Turbo Pascal-Programmen. B.I. Wissenschaftsverlag Mannheim1988, 1991.

[42] Engeln-Mullges, G.; Reutter, F.: Numerik-Algorithmen mit ANSI C-Programmen.(auch fur Turbo Pascal, FORTRAN). B.I. Wissenschaftsverlag Mannheim 1993.

[43] Kose, K.; Schroder, R.; Wieliczek, K.: Numerik sehen und verstehen. Ein kombiniertesLehr- und Arbeitsbuch mit Visualisierungssoftware. Vieweg Braunschweig 1992.

[44] Dietel, J.: Formelsammlung zu Numerischen Mathematik mit Turbo Pascal-Programmen (TPNUM). Rechenzentrum der RWTH Aachen 1993.

[45] Neundorf, W.; Walther, B.: Praktikumsaufgaben fur die Fachgebiete Algorithmen undProgrammierung, Grundlagen der Informatik und Informatik und Programmierung.Skriptum TU Ilmenau 1997.

[46] Schilder, F.: Ein portables C++ Programm zu numerischen Intgeration von Anfangs-wertproblemen. Preprint M 1/97 TUI Marz 1997.

Anschrift:

Dr. Werner NeundorfTechnische Universitat Ilmenau, Institut fur MathematikPF 10 0565D - 98684 Ilmenau

e-mail : [email protected]

59

Gewöhnliche Differentialgleichungen - tu-ilmenau.de · Viele praktische Beispiele führen auf...

Documents

Transcript of Gewöhnliche Differentialgleichungen - tu-ilmenau.de · Viele praktische Beispiele führen auf...