Grafisches Lösen - Ingenieurmathematik (i939)

Um einen ersten Eindruck von der Gestalt der Lösungen einer Differentialgleichung zu bekommen, kann diese grafisch gelöst werden. Dies ist möglich bei

Gleichungen erster Ordnung (2d-Grafik),
Systemen erster Ordnung mit zwei Gleichungen (3d-Grafik),
autonomen Systemen erster Ordnung mit zwei Gleichungen (2d-Grafik),
autonomen Systemen erster Ordnung mit drei Gleichungen (3d-Grafik).

Dabei ist die Arbeit mit 3d-Grafiken allerdings praktisch meist zu aufwendig und zu unübersichtlich.

6.4.1Idee¶

Liegt die zu lösende Differentialgleichung in expliziter Form

y'=f(x,y)

(1)

mit gegebenem $f:\bbR^2\to\bbR$ vor, so kennt man an jedem Punkt $(x,y(x))$ einer Lösung $y$ auch deren Anstieg, denn dieser ist gerade $f(x,y(x))$ . Für jeden Punkt $(x_0,y_0)$ der Ebene kann man also eine Aussage der Form

Wenn eine Lösung durch den Punkt $(x_0,y_0)$ verläuft, so hat sie dort den Anstieg $f(x_0,y_0)$ .

machen. Visualisiert man diesen Anstieg an hinreichend vielen Punkten, so bekommt man einen guten Eindruck vom Verlauf der Lösungskurven.

Die Abbildung

(x,y)\mapsto{\begin{bmatrix}1\\f(x,y)\end{bmatrix}},\quad (x,y)\in\bbR^2

(2)

heißt Richtungsfeld der Differentialgleichung.

6.4.2Umsetzung für Gleichungen¶

Geben uns im Rechteck $[a,b]\times[c,d]$ für geeignet gewählte Stützstellenanzahlen $n_x,n_y\in\bbN$ ein regelmäßiges Gitter von Punkten $(x_k,y,l)$ mit

x_k:=k\,\frac{b-a}{n_x},\quad k=0,\ldots,n_x,\qquad y_l:=l\,\frac{b-a}{n_y},\quad l=0,\ldots,n_y

(3)

vor. An jedem Gitterpunkt wird der Vektor

{\begin{bmatrix}1\\f(x_k,y_l)\end{bmatrix}}

(4)

gezeichnet.

Die Länge der Vektoren ist von untergeordnetem Interesse, da diese aus der Richtung folgt (je steiler der Vektor, das größer die Länge). Wichtig ist nur die Richtung der Vektoren, sodass man auch mit normierten Vektoren arbeiten kann.

In Python kann das Richtungsfeld leicht mit matplotlib.pyplot.quiver visualisiert werden. Die Länge der Vektoren wird so skaliert, dass die Darstellung übersichtlich bleibt.

import matplotlib.pyplot as plt
import numpy as np

# Funktion mit zwei Argumenten (Differentialgleichung)
f = lambda x, y: x * y

# Bildausschnitt und Punktdichte im Gitter
a, b = -3, 3
c, d = -2, 2
nx = 15
ny = 15

# x = Matrix der x-Koordinaten aller Gitterpunkte, y analog
x, y = np.meshgrid(
    np.linspace(a, b, nx),
    np.linspace(c, d, ny)
)

# Richtungen (normiert)
length = np.sqrt(1 + f(x, y) ** 2)  # Länge der Vektoren
dir_x = np.ones_like(x) / length
dir_y = f(x, y) / length

# Lösungskurven als Beispiel
solution_x = np.linspace(a, b, 100)#
solution1_y = np.exp(0.5 * solution_x ** 2)
solution2_y = -0.25 * np.exp(0.5 * solution_x ** 2)

# Plot
fig, ax = plt.subplots()

ax.quiver(x, y, dir_x, dir_y, pivot='middle', angles='xy')

ax.plot(solution_x, solution1_y)
ax.plot(solution_x, solution2_y)

ax.set_aspect('equal', 'box')  # gleiche Einteilung beider Achsen (keine Verzerrung)
ax.set_ylim(c - 0.1 * (d - c), d + 0.1 * (d - c))

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Richtungsfeld zu $y\'=x\\,y$ mit zwei Lösungskurven')

plt.show()

6.4.3Umsetzung für Systeme¶

Für ein explizites System aus zwei Differentialgleichungen

{\begin{align} y_1'&=F_1(x,y_1,y_2),\\ y_2'&=F_2(x,y_1,y_2) \end{align}}

(5)

kann das Richtungsfeld im dreidimensionalen Raum dargestellt werden. Analog zum Vorgehen bei einer einzelnen Differentialgleichung wird hier jedem Punkt $(x,y_1,y_2)\in\bbR^3$ der Vektor

{\begin{bmatrix}1\\F_1(x,y_1,y_2)\\F_2(x,y_1,y_2)\end{bmatrix}}

(6)

zugeordnet. Die Lösungskurven folgen dann wieder der durch die Vektoren gegebenen “Strömung”.

Beispiel (Räuber-Beute-Modell, Lotka 1925, Volterra 1926))

Betrachten ein einfaches abgeschlossenen Ökosystem (z.B. Insel), in dem es Beutetiere gibt, die stets genug Nahrung finden, und Raubtiere, die sich ausschließlich von den Beutetieren ernähren. Zu jeder Zeit $t\geq 0$ bezeichnen wir die Anzahl der Beutetiere mit $B(t)$ und die Anzahl der Raubtiere mit $R(t)$ . Folgende Bedingungen könne für die Beschreibung des Zusammenhangs zwischen den Anzahlen von Raub- und Beutetieren formuliert werden:

Ohne Nahrung sterben die Raubtiere mit der Rate $\alpha_1\in(0,1)$ .
Die Vemehrungsrate der Raubtiere ist proportional zum Nahrungsangebot mit Proportionalitätsfaktor $\beta_1>1$ .
Ohne Angriffe durch Raubtiere vermehren sich die Beutetiere mit der Rate $\alpha_2>1$ .
Die Sterberate der Beutetiere durch Raubtierangriffe ist proportional zur Raubtieranzahl mit Proportionalitätsfaktor $\beta_2\in(0,1)$ .

Diese Bedingungen führen für die Bestimmung der Funktionen $B:[0,\infty)\to[0,\infty)$ und $R:[0,\infty)\to[0,\infty)$ auf das Differentialgleichungssystem

{\begin{align} R'(t)&=R(t)\,\bigl(-\alpha_1+\beta_1\,B(t)\bigr),\\ B'(t)&=B(t)\,\bigl(\alpha_2-\beta_2\,R(t)\bigr). \end{align}}

(7)

Stellen das Richtungsfeld des Lotka-Volterra-Modells grafisch dar, um einen Eindruck vom Verhalten der Lösungskurven zu bekommen. Allerdings wird die dreidimensionale Darstellung meist sehr unübersichtlich.

import numpy as np
import plotly.graph_objects as go

# Differentialgleichung
alpha1 = 1
beta1 = 1
alpha2 = 2/3
beta2 = 4/3
f1 = lambda t, R, B: R * (-alpha1 + beta1 * B)
f2 = lambda t, R, B: B * (alpha2 - beta2 * R)

# Bildausschnitt und Punktdichte im Gitter
t_max = 10
cR, dR = 0, 4
cB, dB = 0, 4
nt = 15
nR = 6
nB = 6

# t = Matrix der t-Koordinaten aller Gitterpunkte, R und B analog
t, R, B = np.meshgrid(
    np.linspace(0, t_max, nt),
    np.linspace(cR, dR, nR),
    np.linspace(cB, dB, nB)
)

# Richtungen (normiert)
length = np.sqrt(1 + f1(t, R, B) ** 2 + f2(t, R, B) ** 2)  # Länge der Vektoren
dir_t = np.ones_like(t) / length
dir_R = f1(t, R, B) / length
dir_B = f2(t, R, B) / length

# Plot
fig = go.Figure()
fig.layout.margin = dict(l=0, r=0, b=0, t=0)  # breiten Rand vermeiden
fig.add_trace(
    go.Cone(
        x=t.reshape(-1), y=R.reshape(-1), z=B.reshape(-1),
        u=dir_t.reshape(-1), v=dir_R.reshape(-1), w=dir_B.reshape(-1),
        colorscale='HSV',
        showscale=False,
        sizeref=0.5
    )
)
fig.update_layout(
    scene_xaxis_title_text='t',
    scene_yaxis_title_text='R(t)',
    scene_zaxis_title_text='B(t)',
    #showlegend=False,
    scene_camera_eye=dict(x=-1.5,y=-1.5,z=1.5)
)
fig.show()

Beim Lotka-Volterra-Modell handelt es sich um ein autonomes System von Differentialgleichungen. Das Richtungsfeld ist bezüglich der Zeit also konstant, sodass die dritte Dimension (Zeitachse) kein zusätzlichen Informationen enthält. Der Verlauf der Lösungskurven kann deutlich übersichtlicher aus dem zweidimensionalen Richtungsfeld zur Zeit $t=0$ (oder jeder anderen Zeit) entnommen werden. Dazu wird die Zeitkoordinate der Richtungen einfach ignoriert.

# Differentialgleichung
alpha1 = 0.75
beta1 = 1.25
alpha2 = 0.6
beta2 = 1.4
f1 = lambda R, B: R * (-alpha1 + beta1 * B)
f2 = lambda R, B: B * (alpha2 - beta2 * R)

# Bildausschnitt und Punktdichte im Gitter
cR, dR = 0, 2
cB, dB = 0, 2
nR = 15
nB = 15

# R = Matrix der R-Koordinaten aller Gitterpunkte, B analog
R, B = np.meshgrid(
    np.linspace(cR, dR, nR),
    np.linspace(cB, dB, nB)
)

# Richtungen (normiert)
length = np.sqrt(1 + f1(R, B) ** 2 + f2(R, B) ** 2)  # Länge der Vektoren
dir_R = f1(R, B) / length
dir_B = f2(R, B) / length

# Plot
fig, ax = plt.subplots()

ax.quiver(R, B, dir_R, dir_B, pivot='middle', angles='xy')

ax.set_aspect('equal', 'box')  # gleiche Einteilung beider Achsen (keine Verzerrung)
ax.set_ylim(cR - 0.1 * (dR - cR), dR + 0.1 * (dR - cR))
ax.set_ylim(cB - 0.1 * (dB - cB), dB + 0.1 * (dB - cB))

ax.set_xlabel('R')
ax.set_ylabel('B')
ax.set_title('Richtungsfeld zum Lotka-Volterra-Modell')

plt.show()

Für autonome Systeme mit drei Gleichungen kann man eine entsprechende dreidimensionale Darstellung erstellen, die aber aufgrund der mangelhaften Übersichtlichkeit praktisch kaum relevant ist.

6.4 Grafisches Lösen

6.4.1Idee¶

6.4.2Umsetzung für Gleichungen¶

6.4.3Umsetzung für Systeme¶