12.9    Lösung von Differenzialgleichungen mit SciPy

Das Federpendel besteht aus einer Schraubenfeder und einer Masse m. Das linke Ende der Feder ist an der Wand befestigt, und an ihrem rechten Ende ist sie mit einer Masse verbunden.

Abbildung 12.12     Feder-Masse-Schwinger

Wenn die Masse um die Strecke x₀ nach rechts ausgerichtet und anschließend losgelassen wird, schwingt sie in Richtung der x-Achse hin und her. Die Frequenz der Schwingung ist von der Masse m und der Federkonstanten c abhängig. Die Dämpfung d bestimmt, wie lange der Schwingvorgang andauert.

Das Pendel wird mit der Kraft

beschleunigt. Die Kraft F_B hat ein negatives Vorzeichen, weil sie in die entgegengesetzte Richtung der Spannkraft F_S wirkt.

Der Einfluss der Dämpfung d ist proportional zur Geschwindigkeit:

Die Spannkraft F_S der Feder wird durch das hookesche Gesetz

bestimmt.

Während des Schwingvorgangs muss die Summe der drei Kräfte 0 sein. Die Pendelschwingungen können also durch die DGL

beschrieben werden. Wenn beide Seiten der DGL durch die Masse m geteilt werden, dann erhält man Ausdrücke für die Abklingkonstante 2δ = d/m und die Kreisfrequenz ω₀ der Pendelschwingungen:

Der Ausdruck vor der 1. Ableitung von x wird als Abklingkonstante bezeichnet:

Die Wurzel aus dem Quotienten der Federkonstante c und der Masse m

ist die Kreisfrequenz der Pendelschwingung.

Für die Anzahl der Schwingungen pro Sekunde gilt:

Daraus kann die Periodendauer einer Schwingung

berechnet werden.

Wenn Sie eine DGL 2. Ordnung numerisch lösen wollen, müssen Sie diese DGL in ein System 1. Ordnung umformen:

SciPy stellt die Funktion

solve_ivp(fun,t_span,y0,method='RK45',t_eval=None, dense_output=False,events=None,vectorized=False,args=None, **options)

für die numerische Lösung von Differenzialgleichungen zur Verfügung.

Als erstes Argument wird der Funktion solve_ivp() das DGL-System fun übergeben. Das zweite Argument t_span steht für den Wertebereich des Zeitintervalls, für das die DGL gelöst werden soll. Das dritte Argument y0 enthält die Anfangswerte. Der Parameter method='RK45' (Standardeinstellung) bedeutet, dass die DGL mit dem expliziten Runge-Kutta-Verfahren der Ordnung 5(4) gelöst wird. Dem Parameter args=(c,m,d) werden die physikalischen Parameter (Federkonstante c, Masse m, Dämpfung d) als Tupel zugewiesen. Die Bedeutung des Parameters dense_output wird in der Quelltextanalyse erklärt. Die anderen Parameter sind optional und werden in unserem Beispiel nicht benötigt.

Listing 12.12 zeigt exemplarisch die Implementierung der SciPy-Funktion solve_ivp(). In den Zeilen 12 bis 15 können Sie die Daten des Feder-Masse-Schwingers ändern.

01  #12_scipy_dgl.py
02  import numpy as np
03  import matplotlib.pyplot as plt
04  from scipy.integrate import solve_ivp
05  #DGL-System
06  def dgl(t,xa,c,m,d):
07      x,v = xa
08      dx_dt = v
09      dv_dt =  -c/m*x - d/m*v
10      return dx_dt,dv_dt
11  #Daten
12  m=1.0  #Masse kg
13  d=0.2  #Dämpfung
14  c=9.87 #Federkonstante N/m
15  x0=5   #Anfangswert der Auslenkung
16  tmax=5
17  #Lösen der DGL
18  t = np.linspace(0,tmax,300)
19  z = solve_ivp(dgl,[0,tmax],[x0,0],args=(c,m,d),dense_output=True)
20  x,v = z.sol(t)
21  #Grafikbereich
22  fig,ax = plt.subplots(2,1,label='Feder-Masse-Schwinger')
23  #Auslenkung
24  ax[0].plot(t,x)
25  ax[0].set(ylabel='x',title='Auslenkung')
26  #Geschwindigkeit
27  ax[1].plot(t,v)
28  ax[1].set(xlabel='Zeit',ylabel='v',title='Geschwindigkeit')
29  fig.tight_layout()
30  plt.show()

Listing 12.12     Lösung einer DGL mit »solve_ivp()«

Ausgabe

Abbildung 12.13     Auslenkung und Geschwindigkeit eines Feder-Masse-Systems

Analyse

In den Zeilen 06 bis 10 steht die Funktionsdefinition dgl(t,xa,...) für die explizite Form des DGL-Systems 1. Ordnung. Die Parameterliste muss immer mit dem Zeitintervall beginnen, in dem die DGL gelöst werden soll. An zweiter Stelle steht der Parameter xa für die Anfangswerte. Bei beiden Parametern muss es sich um eine Sequenz handeln: Liste oder Tupel.

In Zeile 19 löst die SciPy-Funktion solve_ivp(dgl,[0,tmax],[x0,0],...) die DGL. Als erstes Argument muss der Name dgl übergeben werden. Danach erfolgt die Übergabe des Wertebereichs für das Zeitintervall und der Anfangswerte xa als Liste.

Die numerischen Werte der Lösung werden in die Objektvariable z gespeichert. Mit print(z) können Sie sich diese Lösung ausgeben lassen. Sie werden feststellen, dass die Werte der unabhängigen Variablen t in das NumPy-Array t:[...] gespeichert werden. Die Werte der Auslenkung x und der Geschwindigkeit v werden in das NumPy-Array y:[[...]] gespeichert. Der Parameter dense_output ist wichtig. Ihm muss der boolesche Wert True zugewiesen werden, damit die Werte für die Auslenkung und die Geschwindigkeit der Masse aus dem Array y getrennt werden können. Mit print(type(z.t)) und print(type(z.y)) können Sie den Typ von t und y ermitteln. In beiden Fällen erhalten Sie als Ausgabe <class 'numpy.ndarray'>.

In Zeile 20 trennt die Anweisung z.sol(t) die Lösungen und speichert sie in das Tupel x,v.