Python ist eine Skriptsprache mit der besonders einfach wissenschaftliche Anwendungen programmiert werden können. Ein besonders hervorzuhebendes Leistungsmerkmal von Python ist die Erweiterbarkeit durch Module (Programmbibliotheken). Mit dem Modul NumPy können umfangreiche numerischen Berechnungen durchgeführt werden. Mit dem Modul Matplotlib lassen sich aussagekräftige Funktionsplots erstellen. Das Modul SciPy stellt Methoden zur numerischen Berechnung von Differenzialen, Integralen und Differenzialgleichungen zur Verfügung.
Die Beispiele stammen aus dem Buch Python-Kurs für Ingenieure und Naturwissenschaftler (Rheinwerkverlag).
Das Beispiel zeigt den grundsätzlichen Aufbau eines Python-Programms.
Es berechnet die die Stromstärke in einem einfachen Stromkreis.
\[I=\frac{U}{R}\]
#!/usr/bin/env python3
#01_linear1.py
U=230
R=11.8
I=U/R
b="Die Stromstaerke betraegt:"
print(b, I, " A")
Die Spannung u(t), Stromstärke i(t) und die Leistung p(t) eines 50-Hz-Wechselstromnetzes werden in drei Unterdiagrammen dargestellt.
#drei Unterdiagramme
import numpy as np
import matplotlib.pyplot as plt
f=50
URms=230
R=0.001
Xc=10
XL=0
Z= np.sqrt(R**2 + (XL-Xc)**2)
phi=np.arctan((XL-Xc)/R)
I=URms/Z
t = np.linspace(0.0, 20, 1000)
u=np.sqrt(2)*URms*np.sin(2*np.pi*f*t/1000)
i=np.sqrt(2)*I*np.sin(2*np.pi*f*t/1000-phi)
p=u*i
#drei Unterdiagramme erzeugen
plt.subplot(3,1,1)
#Spannung
plt.plot(t, u,color="blue",linewidth=2)
plt.ylabel('u(t)')
plt.grid(True)
#Strom
plt.subplot(3,1,2)
plt.plot(t,i,color="red",linewidth=2)
plt.ylabel('i(t)')
plt.grid(True)
#Leistung
plt.subplot(3,1,3)
plt.plot(t,p,color="green",linewidth=2)
plt.ylabel('p(t)')
plt.xlabel('Zeit in ms')
plt.grid(True)
plt.tight_layout()
plt.savefig("leistung.svg")
plt.show()
Das Beispiel zeigt, wie ein sinusförmiger Strom, der durch eine Spule fließt, differenziert wird.
Für den Spannungsfall an der Spule gilt:
\[u_{L}=L\frac{di}{dt}\]
Die Ableitung wird mit der SciPy-Methode derivative(...) berechnet.
#Spannungsfall an einer Spule
import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative
L=1 #H
f=50 #Frequenz
omega=2*np.pi*f
imax=1. #A
#Stromquelle
def i(t):
return imax*np.sin(omega*t/1e3)
#induktiver Spannungsfall
def u(t):
return L*1e3*derivative(i,t,dx=1e-6)
t = np.linspace(0,20,500)
#zwei Unterdiagramme erzeugen
plt.subplot(2,1,1)
plt.title("eingepraegter Strom")
plt.plot(t, i(t), color="r", lw=2)
plt.ylabel('Strom in A')
plt.grid(True)
#induktiver Spannungsfall
plt.subplot(2,1,2)
plt.title("induktiver Spannungsfall")
plt.plot(t, u(t), color="b",lw=2)
plt.xlabel('t ms')
plt.ylabel('Spannug in V')
plt.grid(True)
plt.tight_layout()
plt.show()
Das Beispiel zeigt, wie die Spannung an einem Kondensator, der an eine Konstantstomquelle angeschlossen ist, berechnet wird.
Für den Spannungsfall am Kondensator gilt:
\[u_{C}=\frac{1}{C} \int^{t}_{0} id\tau \]
Das Integral wir mit der SciPy-Methode quad(...) berechnet.
#Konstante integrieren
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral
C=1 #F
imax=10
#Konstantstromquelle
@np.vectorize
def i(t):
return imax
#Kondensatorspannung
@np.vectorize
def u(t):
uc=(1/C)*integral.quad(i,0,t)[0]
return uc
x = np.linspace(0, 20, 500)
plt.figure(figsize=(8,8))
plt.subplot(2,1,1)
plt.title('Kondensatorstrom')
plt.plot(x, i(x), color="r", lw=2)
plt.ylabel('Strom in A')
plt.grid(True)
#kapazitiver Spannungsfall
plt.subplot(2,1,2)
plt.title('Spannung am Kondensator')
plt.plot(x, u(x), color="b",lw=2)
plt.xlabel('t in s')
plt.ylabel('$u_c(t)$ in V')
plt.grid(True)
plt.show()
Das Beispiel zeigt, wie für ein R-L-Glied die DGL 1.Ordnung gelöst wird. Berechnet wird der Strom iL=f(t) für eine sprungförmige Erregung.
Aus der Schaltung kann folgende DGL abgelesen werden:
\[\frac{di}{dt} =\frac{U_{0}-R\cdot i}{L}\]
Die DGL wird mit der SciPy-Methode solve_ivp(...) gelöst.
#Lineare DGL 1. Ordnung
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
U0=10
R,L=1,1
Tau=L/R
tmax=5*Tau
ti=[0,tmax] #Integrationsintervall
IL0 =[0] #Anfangswert
#DGL 1. Ordnung
def dgl(t,ia):
i=ia #Anfangswert
di_dt=(U0-R*i)/L
return di_dt
t = np.linspace(0,tmax,500)
#Loesung der DGL
z = solve_ivp(dgl,ti,IL0,dense_output=True)
uL_t = z.sol(t)
#Darstellung der Loesung
plt.xlabel("t")
plt.ylabel("$I_L(t)$")
plt.grid(True)
plt.plot(t, uL_t.T,"r",lw=2)
plt.savefig("dgl.svg")
plt.show()
