Differenzialgleichungen

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).

Eine DGL erster Ordnung hat die explizite Form:

\[\ \frac{du}{dt} =f\left( t,u\left( t \right) \right) \]

Eine DGL n-ter Ordnung hat die explizite Form:

\[\ \frac{d^{n}u}{dt^{n}} =f\left( t,u,\frac{du}{dt} ,\ ...\ ,\frac{d^{n-1}u}{dt^{n-1}} \right) \]

Damit Differnzialgleichungen höherer Ordnung mit Python und Julia numerisch gelöst werden können, müssen sie in ein DGL-System erster Ordnung umgewandelt werden.

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.

Symbolsche Lösung mit SymPy: Reihenschwingkreis

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()

Ausgabe

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.

Numerische Lösung: Das Euler-Verfahren

Um die numerische Lösungsmethoden von Python und Julia besser zu verstehen, ist es für den Lernprozess hilfreich sich zunächst das einfachste numerische Lösungsverfahren, das sogennante Euler-Verfahren, genauer anzuschauen. Es arbeitet mit folgendem Algorithmus:


Eingabe a   #untere Integrationsgrenze
Eingabe b   #obere Integrationsgrenze
h = 1e-3    #Zeitschrittweite festlegen
n = (b-a)/h #Anzahl der Rechenschritte berechnen
t0=a
y=y0        #Anfangswert
wiederhole i=1 bis n
    t = t0 + i*h         #Summenalgorithmus
    y = y + h*dgl(t,y)   #Summenalgorithmus und Funktionsaufruf der DGL
gebe t und y als Array zurück

In der Tabelle werden die Implementierungen für Python und Julia gegenüber gestellt. Gelöst wird die DGL y'(t) = y(t). Die Lösung ist die e-Funktion. Denn die Ableitung der e-Funktion ist wieder die e-Funktion. Die numerische Lösung wird mit der analytischen Lösung verglichen.

Python Julia

Python	Julia
#dgl_euler.py import matplotlib.pyplot as plt #Funktionsdefinition für Euler-Verfahen def euler(dgl,y0,a,b,h=1e-3): n=int((b-a)/h) t0=a y=y0 lt,ly=[],[] for i in range(n): t = t0 + ih y = y + hdgl(t,y) lt.append(t) ly.append(y) return lt,ly #Daten a=0 b=2 y0=1 #Anfangswert h=1e-1 #Zeitschrittweite #Definition der DGL def dgl(t,y): return y #1 #Funktionsaufruf vt,vy=euler(dgl,y0,a,b,h) #Grafikbereich fig,ax=plt.subplots() ax.plot(vt,vy,lw=1.5,label="numerisch") ax.grid() ax.set(xlabel="t",ylabel="y") ax.set_title("Euler-Verfahren") #genau import numpy as np t=np.arange(a,b,1e-2) ax.plot(t,np.exp(t),"r--",lw=3, label="genau") ax.legend() plt.show()	#dgl_euler.jl using Plots #Funktionsdefinition für Euler-Verfahen function euler(dgl,y0,a,b,h=1e-3) n = Int64((b - a)/h) t0=a y=y0 vt = Vector{Float64}() vy = Vector{Float64}() for i in 0:n t = t0 + i * h y = y + h * dgl(t,y) push!(vt, t) push!(vy, y) end return vt,vy end #Daten a = 0.0 b = 2.0 y0= 1.0 #Anfangswert h=1e-1 #Zeitschrittweite #Definition der DGL dgl(t,y) = y #1 #Funktionsaufruf vt, vy = euler(dgl,y0,a,b,h) #Grafikbereich plot(vt,vy,xlabel="t",ylabel="y", lw=1.5,label="numerisch") title!("Euler-Verfahren") #genau plot!(a:1e-2:b,t->y0*exp(t), ls=:dash,lw=3,label="analytisch")
Ausgabe	Ausgabe


#dgl_euler.py
import matplotlib.pyplot as plt
#Funktionsdefinition für Euler-Verfahen
def euler(dgl,y0,a,b,h=1e-3):
    n=int((b-a)/h)
    t0=a
    y=y0
    lt,ly=[],[]
    for i in range(n):
        t = t0 + i*h
        y = y + h*dgl(t,y)
        lt.append(t)
        ly.append(y)
    return lt,ly
#Daten
a=0
b=2
y0=1    #Anfangswert
h=1e-1  #Zeitschrittweite
#Definition der DGL
def dgl(t,y):
    return y     #1
#Funktionsaufruf
vt,vy=euler(dgl,y0,a,b,h)
#Grafikbereich
fig,ax=plt.subplots()
ax.plot(vt,vy,lw=1.5,label="numerisch")
ax.grid()
ax.set(xlabel="t",ylabel="y")
ax.set_title("Euler-Verfahren")
#genau
import numpy as np
t=np.arange(a,b,1e-2)
ax.plot(t,np.exp(t),"r--",lw=3,
    label="genau")
ax.legend()
plt.show()


#dgl_euler.jl
using Plots
#Funktionsdefinition für Euler-Verfahen
function euler(dgl,y0,a,b,h=1e-3)
    n = Int64((b - a)/h) 
    t0=a
    y=y0
    vt = Vector{Float64}()
    vy = Vector{Float64}()
    for i in 0:n
        t = t0 + i * h
        y = y + h * dgl(t,y)
        push!(vt, t)
        push!(vy, y)
    end
    return vt,vy
end
#Daten
a = 0.0
b = 2.0
y0= 1.0  #Anfangswert
h=1e-1   #Zeitschrittweite
#Definition der DGL
dgl(t,y) = y    #1
#Funktionsaufruf
vt, vy = euler(dgl,y0,a,b,h)
#Grafikbereich
plot(vt,vy,xlabel="t",ylabel="y",
    lw=1.5,label="numerisch")
title!("Euler-Verfahren")
#genau
plot!(a:1e-2:b,t->y0*exp(t),
    ls=:dash,lw=3,label="analytisch")

Ausgabe Ausgabe

Die Zeitschrittweite wurde bewußt relativ groß gewählt, damit man den Fehler mit bloßem Auge erkennen kann. In #1 können Sie andere DGL-Funktionen eintragen. Für praktische Zwecke eignet sich das Euler-Verfahren nicht. Es leistet aber gute Dienste im Unterricht und in Selbstlernphasen.

Fadenpendel

An einem Faden wird eine Kugel mit der Masse m befestigt. Die Masse wird um den Winkel 𝜑 ausgelenkt und losgelassen. Das Pendel beginnt zu schwingen. Übrigens: Das foucaultsches Pendel hat eine Länge von 67 Meter und eine Masse von 28 Kilogramm.

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.

Lösung mit Python

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()

Ausgabe

Analyse

In #1 wird das DGL-System definiert. In #2 wird es gelöst und in #3 werden die Lösungen getrennt. Die Standardintegrationsmethode ist das Runge-Kutta-Verfahren RK45. Es können der SciPy-Funktion solve_ivp() auch andere Integrationsmethoden wie z.B.: RK23, DOP853, Radau, BDF und LSODA. Die absolute und die relative Genauigkeit können mit atol und rtol vorgegeben werden.

Lösung mit Julia

Julia löst Differenzialgleichungen mit den Funktionen ODEProblem(...) und solve(...).

Falls noch nicht erledigt, müssen Sie das Paket DifferentialEquations installieren.


using Pkg 

Pkg.add("DifferentialEquations")

Die Pendellänge und die Dämpfungskonstante sind als globale Variablen deklariert. Sie können auch als Parameter p = [l,d] in einen Vektor gespeichert werden und in #2 der Funktion ODEProblem(fadenpendel, u0, tspan,p) als Parameter übergeben werden.


#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"])

Ausgabe

Analyse

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. Die Funktion ODEProblem() beschreibt "was" gelöst werden soll, während die Funktion solve() bestimmt, "wie" es gelöst werden soll. Der solve-Funktion können noch Integrationsmethoden, wie z.B.: Midpoint(), Heun(), RK4(), Tsit5() ... als Argumente übergeben werden. Die Standardintegrationsmethode ist eine Hermite-Polynominterpolation 3. Ordnung. Außerdem kann man auch die relative Genauigkeit reltol und die absolute Genauigkeit abstol vorgeben.

Sprungantwort eines fremderregten Gleichstrommotors

Das Beispiel zeigt, wie die DGL für die Sprungantwort eins R-L-C-Gliedes 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 = C_dyn wird aus dem Trägheitsmoment J, dem Bemessungsstrom I_n und dem Bemessungsmoment M_n 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.

Lösung mit Python

Die DGL wird mit der SciPy-Funktion solve_ivp(...) gelöst.

#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()

Ausgabe

sprunganwort gleichstrommotor

Analyse

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.

Lösung mit Julia

Die DGL wird mit den Funktionen ODEProblem(...) und solve(...) gelöst.


#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)

Ausgabe

sprunganwort gleichstrommotor

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. Dieses DGL-System ist nicht-linear, es kann nur numerisch gelöst werden.

Lösung mit Python

#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()

Ausgabe

Analyse

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.

Lösung mit Julia


#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"])

Ausgabe

Besuchen Sie die Seite von SciPy.
Besuchen Sie die Seite von Julia.