Differenzialgleichungen mit SciPy lösen

Differenzialgleichungen werden mit der Methode
z=solve_ivp(dgl,intervall,anfangswerte,args=(...),dense_output=True)
gelöst. Der Lösungsvektor wird in dem Objekt z gespeichert. Die einzelnen Lösungen erhält man durch slicing. Differenzialgleichungen n-ter Ordnung müssen in ein DGL-System 1. Ordnung umgewandelt werden.

Sprungantwort eines fremderregten Gleichstrommotors

Das Beispiel zeigt, wie die DGL für die Sprungantwort eins R-L-C-Gliedes mit der SciPy-Methode solve_ivp(...) numerisch gelöst wird. Das R-L-C-Glied kann als Ersatzschaltbild eines fremderregten Gleichstrommotors interpretiert werden.

ersatzschaltbild gleichstrommotor

Aus dem Ersatzschaltbild des Gleichstrommotors kann folgendes DGL-System abgelesen werden:

\[\frac{du_{c}}{dt} =\frac{i}{C} \] \[\frac{di}{dt} =\left( U_{0}-RC\frac{du}{dt} -u_{c}\right) \frac{1}{L} \]

Die dynamische Kapazität C = Cdyn wird aus dem Trägheitsmoment J, dem Bemessungsstrom In und dem Bemessungsmoment Mn des Gleichstrommotors berechnet:

\[C=J\left( \frac{I_{n}}{M_{n}} \right)^{2}\]

Sie verzögert die Hochlaufzeit.

Die elektrische Zeitkonstante
\[T_{el}=\frac{L}{R} \] beeinflusst die Verzögerung des Stromflusses während des Einschaltvorgangs.

Aus dem Ankerwiderstand R und der dynamischen Kapazität C wird die mechanische Zeitkonstante
\[T_{mech}=R\cdot C\] berechnet.

#Lineare DGL 2. Ordnung
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
U0 = 100  #Eingangsspannung
R = 1.5   #Ankerwiderstand
L = 0.025 #Ankerinduktivität
Mn=150    #Bemessungsmoment
In=50     #Bemessungsstrom
J=0.2     #Trägheitsmoment
tmax=0.5
#DGL-System
def dgl(t,anfangswerte,R,L,C):
    uc,i = anfangswerte
    duc_dt = i/C
    di_dt = (U0 - R*C*duc_dt-uc)/L
    return [duc_dt, di_dt]

C = J*(In/Mn)**2 #dynamische Kapazität
a0 = [0,0]       #Anfangswerte
ti=[0,tmax]      #Integrationsintervall
t = np.linspace(0,tmax,500)
#Lösung des DGL-Systems 
z=solve_ivp(dgl,ti,a0,args=(R,L,C),dense_output=True) #1
uc_i_t = z.sol(t)
#Trennen der Lösung
uc_t = uc_i_t.T[ : , 0]
i_t  = uc_i_t.T[ : , 1]
#zwei Unterdiagramme 
fig,axes=plt.subplots(2,1,figsize=(8,8))
#Kondensatorspannung
axes[0].set_title("Sprungantwort eines fremderregten Gleichstrommotors")
axes[0].plot(t, uc_t,"b",lw=2)
axes[0].set_ylabel('Ausgangsspannung in V')
axes[0].grid(True)
#Stromverlauf
axes[1].plot(t, i_t,"r",lw=2)
axes[1].grid(True)
axes[1].set_ylabel('Ankerstrom in A')
axes[1].set_xlabel('Zeit in Sekunden')
print("dynamische Kapazität:",C,"F")
fig.tight_layout()
plt.show()

Ausgabe

dynamische Kapazität: 0.022222222222222223 F

sprunganwort gleichstrommotor

Analyse

Die Lösung der DGL erfolgt in Zeile #1 mit der SciPy-Methode solve_ivp(dgl,ti,a0,args=(R,L,C),dense_output=True). Als erster Parameter wird die Funktionsdefinition der DGL übergeben. Der Parameter ti legt das Integrationsintervall fest. Der dritte Parameter a0 bestimmt die Anfangswerte. Mit args=(R,L,C) können zusätzliche Parameter übergeben werden. Der letzte Parameter dense_output=True ermöglicht, dass die Lösung der DGL in dem Objekt z gespeichert werden kann.

Simulation einer Epidemie

Die Ausbreitung einer Epidemie lässt sich mit einem Differenzialgleichungssystem aus drei Differenzialgleichungen beschreiben. Eine Population N wird in drei Gruppen aufgeteilt: Gesunde S, Infizierte I und wieder Genesene R.

\[\frac{dS}{dt} =-b\frac{S\cdot I}{N} \]

\[\frac{dI}{dt} =b\frac{S\cdot I}{N} -g\cdot I\]

\[\frac{dR}{dt} =g\cdot I\]

Dabei steht g für die Genesungsrate und b für die Infektionsrate.

#Simulation einer Epidemie
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
tmax=120
#Anfangswerte
S0=997     #nicht immune Gesunde
I0=3       #Infizierte
R0=0       #Genesene
N=S0+I0+R0 #Population
#Parameter
b=0.4  #Infektionsrate
g=0.04 #Genesungsrate
#DGL-System
def dgl(t,ya):
    S,I,R=ya
    dS_dt=-b*S*I/N    #nicht immune Gesunde
    dI_dt=b*S*I/N-g*I #Infizierte
    dR_dt=g*I         #Genesene
    return [dS_dt,dI_dt,dR_dt]
#Anfangswerte
y0 = [S0,I0,R0]
ti=[0,tmax]
t = np.linspace(0, tmax, 500)
z=solve_ivp(dgl,ti,y0,dense_output=True)
y_t = z.sol(t)
#Trennen der Lösung
S = y_t.T[:,0] #1
I = y_t.T[:,1] #2
R = y_t.T[:,2] #3
plt.figure(figsize=(8,6))
plt.plot(t, S,'b',lw=2,label="Gesunde")
plt.plot(t, I,'r',lw=2,ls="dashed",label="Infizierte")
plt.plot(t, R,'g',lw=2,ls="dotted",label="Genesene")
plt.legend(loc="best")
plt.xlabel("Zeit")
plt.ylabel("Individuen")
plt.grid(True)
plt.show()        

Ausgabe

epidemie simulation

Analyse

In den Zeilen #1 bis #3 werden durch slicing und Transponierung T die einzelnen Lösungsvektoren für die Kranken (S), die Infizierten (I) und die wieder Genesenen (R) ermittelt.

Besuchen Sie die Seite von SciPy.

Startseite NumPy Matplotlib SymPy VPython