SciPy stellt Funktionen für das numerische Differenzieren und Integrieren bereit.
Außerdem können Sie mit diesem Modul Differenzialgleichungen und Differenzialgleichungssysteme numerisch lösen.
Mit den SciPy-Funktionen derivative() können mathematische Funktionen numerisch differenziert werden.
Das Untermodul misc wird demnächst aus dem Modul scipy entfernt. Damit steht auch die Funktion derivative() nicht mehr zur Verfügung. Stattdessen empfiehlt die Dokumentation von SciPy das Modul numdifftools zu verwenden.
#ableitung.py
import numpy as np
import matplotlib.pyplot as plt
from numdifftools import Derivative
#Funktionsdefinition
def f(x):
#y=5*x
y=0.25*x**2
#y=np.sin(x)
return y
def ableitung(x,n=1):
df = Derivative(f,n=n)
return df(x)
n=1 #Ordnung der Ableitung
x=np.linspace(0,10,100)
fig, ax=plt.subplots(figsize=(5,5))
ax.plot(x,f(x),'b-',label='Funktion')
ax.plot(x,ableitung(x,n),'r--',label='Ableitung')
ax.legend(loc='best')
ax.set(xlabel='x',ylabel='y',title='Ableitung einer Funktion')
fig.savefig('ableitung.svg')
plt.show()
Mit den SciPy-Funktionen quad(), dblquad() und tplquad() können Sie Einfachintegrale, Zweifachintegrale und Dreifachintegrale numerisch berechnen.
Das Programm berechnet die Stammfunktion einer stetigen mathematischen Funktionen.
\[\frac{1}{4} \int x^{2}dx=\frac{1}{12} x^{3} \]
#integral.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
#Funktionsdefinition
def f(x):
#y=5*x
y=x**2/4
#y=np.sin(x)
return y
@np.vectorize
def integral(x):
return quad(f,0,x)[0]
x=np.linspace(0,10,100)
fig, ax=plt.subplots(figsize=(5,5))
ax.plot(x,f(x),'b-',label='Funktion')
ax.plot(x,integral(x),'r--',label='Integral')
ax.legend(loc='best')
ax.set(xlabel='x',ylabel='y',title='Ableitung einer Funktion')
fig.savefig('integral.svg')
plt.show()
\[\int^{2}_{1} \int^{3y}_{y} \left( x+y\right) dxdy=14 \]
#zweifach_dblquad.py
from scipy.integrate import dblquad
#Funktionsdefinition
def f(x,y):
return x+y
#untere Grenze, variabel
def y1(y):
return y
#obere Grenze, variabel
def y2(y):
return 3*y
#Hauptprogramm
V= dblquad(f,1,2,y1,y2)[0]
print("V =",V)
V = 14.0
\[\int^{1}_{0} \int^{1-x}_{0} \int^{2-x}_{0} xyzdzdydx=\frac{13}{240} \]
#deifachintegral.py
from scipy.integrate import tplquad
#Funktionsdefinition
def f(z,y,x):
return x*y*z
#obere Grenze, variabel
def y1(x):
return 1-x
#obere Grenze, variabel
def y2(x,y):
return 2-x
#Hauptprogramm
I = tplquad(f,0,1,0,y1,0,y2)[0]
print("I =",I)
I = 0.054166666666666675
Differenzialgleichungen werden mit der Funktion
z=solve_ivp(dgl,intervall,anfangswerte,args=(...),dense_output=True)
numerisch gelöst. Der Lösungsvektor wird in dem Objekt z gespeichert. Die einzelnen Lösungen werden mit x, y = z.sol(t) in das Tupel x,y gespeichert. Differenzialgleichungen n-ter Ordnung müssen in ein DGL-System 1. Ordnung umgewandelt werden.
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.
#Lineare DGL 2. Ordnung
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
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]
#Berechnungen und Anfangswerte
C = J*(In/Mn)**2 #dynamische Kapazität
a0 = [0,0] #Anfangswerte fuer Kondensatorspannung und Strom
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,ic = z.sol(t)
#Grafikbereich
fig,axes=plt.subplots(2,1,figsize=(8,8))
#Kondensatorspannung
axes[0].set_title("Sprungantwort eines fremderregten Gleichstrommotors")
axes[0].plot(t, uc,"b",lw=2)
axes[0].set_ylabel('Ausgangsspannung in V')
axes[0].grid(True)
#Stromverlauf
axes[1].plot(t, ic,"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()
dynamische Kapazität: 0.022222222222222223 F
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.
Die Massen m von zwei Fadenpendeln werden mit einer Schraubenfeder verbunden. Die Schwingungen sind von den Massen, der Federkonstanten c und den Fadenlängen l abhängig.
Das gekoppelte Fadenpendel kann durch folgendes Differenzialgleichungssystem beschrieben werden:
\[\frac{d\varphi_{1} }{dt} =\omega_{1}\]
\[\frac{d\omega_{1} }{dt} =-\frac{g}{l} \varphi_{1} -\frac{c}{m} \left( \varphi_{1} -\varphi_{2} \right) \]
\[\frac{d\varphi_{2} }{dt} =\omega_{2}\]
\[\frac{d\omega_{2} }{dt} =-\frac{g}{l} \varphi_{2} +\frac{c}{m} \left( \varphi_{1} -\varphi_{2} \right) \]
Das DGL-System kann direkt in eine Python-Funktion implementiert werden. Gelöst wird das DGL-System wieder mit der SciPy-Methode solve_ivp(...).
#gekoppeltes_fadenpendel.py
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
#Pendeldaten
g = 9.81 #Erdbeschleunigung in m/s^2
l= 0.2 #Pendellaengen in m
m= 0.2 #Pendelmassen in kg
c=1 #Federkonstante N/m
phi1, phi2 = -20, 0 #1
tmax=10
#DGL-System, der Unterstrich bedeutet 1. Ableitung nach der Zeit
def dgl(t,ya,l, m, c):
phi1,w1,phi2,w2 = ya #1
phi1_dt=w1
w1_dt=-g/l*phi1-c/m*(phi1-phi2)
phi2_dt=w2
w2_dt=-g/l*phi2+c/m*(phi1-phi2)
return [phi1_dt, w1_dt, phi2_dt, w2_dt] #2
#Lösung des DGL-Systems
omega1 = omega2 = 0
ya =[np.radians(phi1),omega1,np.radians(phi2),omega2] #Anfangswerte
t = np.linspace(0,tmax,500)
z=solve_ivp(dgl,[0,tmax],ya,args=(l,m,c),dense_output=True)#3
phi1,w1,phi2,w2 =z.sol(t) #4
auslenkung1 = np.degrees(phi1) #5
auslenkung2 = np.degrees(phi2) #6
#Grafik-Bereich
fig,ax=plt.subplots(2,1,figsize=(6,8))
#Auslenkung l1
ax[0].plot(t,auslenkung1,"r",lw=2)
ax[0].set(ylabel=r"$\varphi_{1}$",title="Pendel 1")
#Auslenkung l2
ax[1].plot(t,auslenkung2,'b',lw=2)
ax[1].set(xlabel="Zeit",ylabel=r"$\varphi_{2}$",title="Pendel 2")
ax[0].grid(True);ax[1].grid(True)
fig.tight_layout()
plt.show()
Zeile #1 legt die Anfangswerte für die Auslenkwinkel fest.
Die Python-Funktion gibt in Zeile #2 das DGL-System als Liste zurück.
In Zeile #3 wird das DGL-System gelöst.
In Zeile #4 wird die Lösung getrennt. Die einzelnen Teillösungen werden in ein Tupel gespeichert.
In den Zeilen #5 und #6 werden die Auslenkwinkel von Radiant in Grad umgerechnet.
Dieses Beispiel wird in meinem Buch "Mathematische Algorithmen mit Python" mit dem Euler-Verfahren gelöst.
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.
#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]
t = np.linspace(0, tmax, 500)
z=solve_ivp(dgl,[0,tmax],y0,dense_output=True)
S, I, R = z.sol(t) #1
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(t, S,'b-', lw=2, label="Gesunde")
ax.plot(t, I,'r--',lw=2, label="Infizierte")
ax.plot(t, R,'g:', lw=2, label="Genesene")
ax.legend(loc="best")
ax.set_xlabel("Zeit")
ax.set_ylabel("Individuen")
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.