Das deutsche Python-Forum

Code: Alles auswählen

import numpy as np
from scipy import integrate
from scipy.integrate import ode
import sympy
import matplotlib.pylab as plt
import math

# Definiere SymPy symbols für die Variablen und die Funktionen
t, l_1, l_2, l_12, c_1, r_2 = sympy.symbols("t, l_1, l_2, l_12, c_1, r_2")   
q1, q2 = sympy.symbols("q1, q2", cls=sympy.Function)

#ODE´s für primärer und sekundärer Schwingkreis
ode1 = sympy.Eq(l_1 * q1(t).diff(t,t) + l_12 * q2(t).diff(t,t) + q1(t)/c_1) 
ode2 = sympy.Eq(l_12 * q1(t).diff(t,t) + l_2 *q2(t).diff(t,t) + r_2 * q2(t).diff(t)) 

# Jetzt sind ode1 und ode2 sympy expressions für die Differentialgleichungen 2. Ordnung (noch nicht in Form für ODE solvers)

# Umschreiben in ein System von DGL´s 1. Ordnung
# z_1 = q1(t); z_2 = q1'(t); z_3 = q2(t), z_4 = q2'(t)

z_1, z_2, z_3, z_4 = sympy.symbols("z_1, z_2, z_3, z_4", cls=sympy.Function)
varchange = {q1(t).diff(t,t):
            z_2(t).diff(t),
            q1(t): z_1(t),
            q2(t).diff(t,t):
            z_4(t).diff(t),
            q2(t): z_3(t)}

ode1_vc = ode1.subs(varchange)
ode2_vc = ode2.subs(varchange)

# wir müssen noch zwei weitere ODE's einführen: z_1'(t) und z_3'(t)
ode3 = z_1(t).diff(t) - z_2(t)
ode4 = z_3(t).diff(t) - z_4(t)
# An dieser Stelle haben wir 4 gekoppelte DGL's 1.Ordnung

# Nun müssen wir die Ableitungen dieser Funktionen lösen um die ODE's in Standardform zu erhalten (benutze sympy.solve)
z = sympy.Matrix([z_1(t), z_2(t), z_3(t), z_4(t)])
vcsol = sympy.solve((ode1_vc, ode2_vc, ode3, ode4), z.diff(t), dict=True)
f = z.diff(t).subs(vcsol[0])
# Nun ist f sympy expression für die ODE Funktion f(t,z(t))
# Wir könnten die ODE mit sympy.Eq(z.diff(t), f) darstellen, aber das Ergebnis wäre sehr lange
# Main purpose: Konstruiere f an dieser Stelle um es in eine NumPy-aware Funktion zu konvertieren, welche mit integrate.odeint oder integrate.ode genutzt werden kann
# Konstruiere so eine Funktion mit sympy.lambdify

k = 0.5 # Kopplungsfaktor \in (0, 1) Induktionen der Spulen müssen sehr weit voneinander entfernt sein => l1 ist VIEL kleiner als l2
params = {l_1: 5, l_2: 100, l_12: 12, c_1: 5, r_2: 10} # Werte der Parameter Anmerkung: Integrieren fehlt!
_f_np = sympy.lambdify((t, z), f.subs(params), 'numpy')
f_np = lambda _t, _z, *args: _f_np(_t, _z)

# benutze jac für Jacobianmatrix
jac = sympy.Matrix([[fj.diff(zi) for zi in z] for fj in f])
_jac_np = sympy.lambdify((t, z), jac.subs(params), 'numpy')
jac_np = lambda _t, _z, *args: _jac_np(_t, _z)
# an dieser Stelle haben wir spezielle Werte der Systemparameter substituiert bevor wir sympy.lambdify aufrufen. 

# vorerst u_omega = 2 !INTEGRIEREN!
z0 = [c_1*2, 0, 0, 0]
tt = np.linspace(0, 20, 1000)
r = integrate.ode(f_np,jac_np).set_initial_value(z0, tt[0])
dt = tt[1] - tt[2]
zz = np.zeros((len(tt)(z0)))
idx = 0
while r.successful() and r.t < tt[-1]:
    zz[idx, :] = r.z
    r.integrate(r.t + dt)
    idx += 1
# die Lösung der ODE's ist nun in dem Array zz gepeichert, welches die Größe (1000, 4) hat. 

# plot in the x-y plane
q1_np, q2_np = zz[:,0], zz[:,2]
x1 = params[l_1] * q1_np
x2 = params[l_2] * q2_np

fig = plt.figure(figsize=(10,4))
ax1 = plt.subplot2grid((2,5),(0,0),colspan=3)
ax2 = plt.subplot2grid((2,5),(1,0),colspan=3)

ax1.plot(tt, x1, 'r')
ax1.set_ylabel('$x_1$')
ax1.set_yticks([-3, 0, 3])

ax2.plot(tt, x2, 'r')
ax2.set_ylabel('$x_2$', fontsize=18)
ax2.set_xlabel('$t$', fontsize=18)
ax2.set_yticks([-3, 0, 3])