Differenzialgleichungen sind grundlegend für die mathematische Modellbildung physikalischer, technischer und sozialer Systeme. Sie beschreiben deren Verhalten und erlauben präzise Vorraussagen ihres Verhaltens. Technische System können durch Simulationen und Variationen einzelner Parameter optimiert werden. In der Praxis erfolgt ihre Lösung in der Regel mit numerischen Methoden (z.B. das Runge-Kutta-Verfahren). Am Beispiel eines Reihenschwingkreises, eines Fadenpendels, dem Einschalten eines Gleichstrommotors und der Simulation einer Epidemie soll gezeigt werden, wie mit Python und Julia Differenzialgleichungen symbolisch und numerisch gelöst werden können.
Die Reihenschaltung aus Induktivität, ohmscher Widerstand und Kapazität wird an ein Stromquelle angeschlossen. Zum Zeitpunkt t = 0 wird der Schalter geschlossen. Untersucht werden soll der Strom i = f (t).
Für die Summe aller Spannungen gilt nach der Maschenstromregel:
\[\ u_{L}+u_{R}+u_{c}=U \]
Berücksichtig man den induktiven und den kapazitiven Spannungsfall erhält man:
\[\ L\frac{di}{dt} +R\cdot i+\frac{1}{C} \int idt=U \]
Durch Ableitung erhält man eine DGL zweiter Ordnung:
\[\ L\frac{d^{2}i}{dt^{2}} + R\frac{di}{dt} +\frac{1}{C} i=0 \]
Beide Seiten der DGL werden mit C multipliziert:
\[\ LC\frac{d^{2}i}{dt^{2}} + RC\frac{di}{dt} + i=0 \]
Die SymPy-Methode dsolve(...) kann nur lineare Differenzialgleichungen (DGLs) symbolisch lösen. Nichtlineare DGLs müssen numerisch gelöst werden. Mit dem SymPy-Programm sym_dgl.py können Sie verschiedene Lösungen testen, indem Sie für die Störfunktion (#3), für die Induktivität L, für die Kapazität C und für den ohmschen Widerstand R andere Werte eintragen.
#sym_dgl.py
from sympy import *
t=symbols('t') #1
i=Function('i')(t) #2
I0=1 #3
L=1
C=1
R=0
dgl=Eq(L*C*i.diff(t,2) + R*C*i.diff(t,1) + i,I0) #4
#Anfangswerte
aw={
i.subs(t,0):0,
i.diff(t,1).subs(t,0):0
}
ia_t=dsolve(dgl,i) #5
is_t=dsolve(dgl,i,ics=aw) #6
i_t=is_t.rhs
#Ausgabe
print("allgemeine Lösung\n",ia_t)
print("spezielle Lösung\n",is_t)
print("rechte Seite der Gleichung\n i(t) =",i_t)
plt=plot(i_t,(t,0,10),show=False,line_color='red',ylabel='i(t)')
#plt.save("plot_rlc_sprungantwort.svg")
plt.show()
allgemeine Lösung Eq(i(t), C1*sin(t) + C2*cos(t) + 1) spezielle Lösung Eq(i(t), 1 - cos(t)) rechte Seite der Gleichung i(t) = 1 - cos(t)
In #1 wird die symbolische Variable t vereinbart. In #2 wird festgelegt, das i eine Funktion von t sein soll. In #3 steht die Störfunktion. In #4 wird die DGL zweiter Ordnung definiert und in #5 wird sie allgemein gelöst. Der Parameter ics=aw (initial/boundary conditions) bewirkt, dass dsolve(...) die spezielle Lösung berechnet.
Aus dem Bild kann folgendes DGL-System abgelesen werden:
\[\frac{d\varphi }{dt} =\omega \]
\[\frac{d\omega }{dt} =-\frac{g}{l} \sin \varphi \]
In der zweiten Gleichung kann noch ein Term für die Dämfung eingefügt werden, sie ist proportional zur Winkelgeschwindigkeit.
Differenzialgleichungen werden mit der SciPy-Funktion
sol=solve_ivp(dgl,intervall,anfangswerte,args=(...), t_eval=t)
numerisch gelöst. Der Lösungsvektor wird in dem Objekt sol gespeichert.
Die einzelnen Lösungen werden mit y1, y2 = sol.y in das Tupel y1,y2 gespeichert.
Differenzialgleichungen n-ter Ordnung müssen in ein DGL-System 1. Ordnung umgewandelt werden.
#dgl_fadenpendel.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
phi0 =1 #Auslenkung in rad
l=1 #Pendellänge in m
d=0.2 #Dämpfungskonstante
a0 = [phi0,0] #Anfangswerte
g=9.81
tmax=5 #Zeit in Sekunden
#DGL-System
def fadenpendel(t,anfangswerte): #1
phi,w = anfangswerte
dphi_dt = w
dw_dt = -d*w -g/l*np.sin(phi)
return [dphi_dt, dw_dt]
#Lösung der DGL
t = np.linspace(0,tmax,500)
sol=solve_ivp(fadenpendel,[0,tmax],a0,t_eval=t) #2
phi,w = sol.y #3
#Grafikbereich
fig,ax=plt.subplots(2,1,figsize=(9,6)) #Objekte erzeugen
ax[0].plot(t,phi,'b-',lw=1.5,label=r'$\varphi$')
ax[1].plot(t,w,'r-',lw=1.5,label=r'$\omega$')
ax[1].set_xlabel('t')
ax[0].legend(loc="best"),ax[1].legend(loc="best")
ax[0].grid(True),ax[1].grid(True)
#fig.savefig('py_fadenpendel.svg')
plt.show()
In #1 wird das DGL-System definiert. In #2 wird es gelöst und in #3 werden die Lösungen getrennt.
Julia löst Differenzialgleichungen mit den Funktionen ODEProblem(...) und solve(...).
using Pkg
Pkg.add("DifferentialEquations")
#dgl_fadenpendel.jl
using Plots, DifferentialEquations
# Daten
phi = 1 #Auslenkung in rad
l = 1.0 #Pendellänge in m
d = 0.2 #Dämpfungskonstante
tspan = (0, 5)
g = 9.81
u0 = [phi, 0.0]
# Definition der DGL
function fadenpendel(du, u, p, t) #1
du[1] = dφ = u[2] #u[2] = ω
du[2] = dω = -d * u[2] - g / l * sin(u[1]) #u[1] = φ
end
# Lösung der DGL
problem = ODEProblem(fadenpendel, u0, tspan) #2
sol = solve(problem) #3
plot(sol, lw=1.5, label=["φ" "ω"], layout=(2, 1), color=[ "blue" "red"])
In #1 wird das DGL-System definiert. Das ODEProblem wird in #2 implementiert, indem zuerst die ODEFunction ODEProblem(fadenpendel, u0, tspan) aufgerufen wird. Die rechte Seite der DGL wird an den Konstruktor übergeben. Das Resultat wird in das Objekt problem gespeichert. In #3 wird die DGL-gelöst. Sie wird in das Objekt sol gespeichert.
Das Beispiel zeigt, wie die DGL für die Sprungantwort eins R-L-C-Gliedes mit der SciPy-Funktion solve_ivp(...)
numerisch gelöst wird. Das R-L-C-Glied kann als Ersatzschaltbild eines fremderregten Gleichstrommotors interpretiert werden.
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_c}{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.
#dgl_gleichstrommotor.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
U0 = 100 #Eingangsspannung in V
R = 1.5 #Ankerwiderstand in Ohm
L = 0.025 #Ankerinduktivität in H
Mn=150 #Bemessungsmoment in Nm
In=50 #Bemessungsstrom in A
J=0.2 #Trägheitsmoment in kgm^2
C = J*(In/Mn)**2 #dynamische Kapazität
a0 = [0,0] #Anfangswerte
tmax=0.5 #Zeit in Sekunden
ti=[0,tmax] #Integrationsintervall
#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]
#Lösung der DGL
t = np.linspace(0,tmax,500)
sol=solve_ivp(dgl,ti,a0,args=(R,L,C),t_eval=t)
uc,ic = sol.y
#Grafikbereich
fig,ax=plt.subplots(2,1,figsize=(9,6)) #Objekte erzeugen
ax[0].plot(t,uc,'b-',lw=1.5,label='Ausgangsspannung')
ax[1].plot(t,ic,'r-',lw=1.5,label='Stromstärke')
ax[1].set_xlabel('t')
ax[0].legend(),ax[1].legend()
ax[0].grid(True),ax[1].grid(True)
plt.show()
Die Lösung der DGL erfolgt in Zeile #1 mit der SciPy-Funktion solve_ivp(dgl,ti,a0,args=(R,L,C),t_eval=t).
Als erster Parameter wird die Funktionsdefinition der DGL übergeben.
Der Parameter ti legt das Integrationsintervall fest.
Der dritte Parameter a0 bestimmt die Anfangswerte für die Stromstärke und Kondensatorspannung. Mit args=(R,L,C) können zusätzliche Parameter übergeben werden.
#dgl_gleichstrommotor.jl
using Plots, OrdinaryDiffEq
U0 = 100 #Eingangsspannung in V
R = 1.5 #Ankerwiderstand in Ohm
L = 0.025 #Ankerinduktivität in H
Mn = 150 #Bemessungsmoment in Nm
In = 50 #Bemsseungsstrom in A
J = 0.2 #Trägheitsmoment
C = J * (In / Mn)^2 #dynamische Kapazität
i0 = [0, 0]
tspan = (0, 0.5)
#DGL für Sprungantwort
function gleichstrommotor(du, u, p, t)
du[1] = duc = u[2] / C #uc(t)
du[2] = di = (U0 - R * C * du[1] - u[1]) / L #ic(t)
end
#ODE Problem
prob = ODEProblem(gleichstrommotor, i0, tspan)
#Lösung
sol = solve(prob)
#Darstellung der Lösung
plot(sol, label=["u(t)" " i(t)"], layout=(2, 1), color=["blue" "red"], lw=1.5)
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. Dieses DGL-System ist nicht-linear, es kann nur numerisch gelöst werden.
#dgl_epidemie.py
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]
t = np.linspace(0, tmax, 500)
sol=solve_ivp(dgl,[0,tmax],y0,t_eval=t)
S, I, R = sol.y #1
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(t, S,'b-', lw=1.5, label="Gesunde")
ax.plot(t, I,'r-',lw=1.5, label="Infizierte")
ax.plot(t, R,'g-', lw=1.5, label="Genesene")
ax.legend(loc="best")
ax.set_xlabel("t")
ax.set_ylabel("N")
ax.grid(True)
plt.show()
In Zeile #1 werden die Teillösungen des Lösungsvektors z in ein Tupel gespeichert.
Das Tupel enthält die Werte für die Gesunden (S: susceptible, ansteckbare), die Infizierten (I: infected,
infizierte) und die wieder Genesenen (R: removed, recovered) als eindimensionales Array.
#dgl_epidemie.jl
using Plots,DifferentialEquations
#Daten
S0 = 997 #nicht imune Gesunde
I0 = 3 #Infizierte
R0 = 0 #Genesene
N = S0 + I0 + R0
u0=[S0,I0,R0]
b=0.4
g=0.04
tspan=(0,120)
# Definition der DGL
function dgl(du,u,p,t)
du[1] = -b*u[1]*u[2]/N #nicht immune Gesunde
du[2] = b*u[1]*u[2]/N - g*u[2] #Infizierte
du[3] = g*u[2] #Genesene
end
# Lösung der DGL
problem = ODEProblem(dgl,u0,tspan)
sol = solve(problem)
plot(sol,lw=1.5,ylabel="N",
label=["Gesunde" "Infizierte" "Genesene"],
color=["blue" "red" "green"])
Besuchen Sie die Seite von SciPy.
Besuchen Sie die Seite von Julia.
Startseite | Funktionsplots | Matrizen | Differenzieren | Integrieren | Animationen |