leider habe ich bei der Verwendung der ODE-Klasse ein Problem.
Folgende Fehlermeldung
PS C:\Users\Julia> & C:/Users/Julia/AppData/Local/Programs/Python/Python36/python.exe c:/Users/Julia/odeGekoppelt.py
File "c:/Users/Julia/odeGekoppelt.py", line 60, in <module>
r = integrate.ode(f_np,jac_np).set_initial_value(z0, tt[0])
File "C:\Users\Julia\AppData\Local\Programs\Python\Python36\lib\site-packages\scipy\integrate\_ode.py", line 366, in set_initial_value
self._y = asarray(y, self._integrator.scalar)
File "C:\Users\Julia\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\core\_asarray.py", line 85, in asarray
return array(a, dtype, copy=False, order=order)
File "C:\Users\Julia\AppData\Local\Programs\Python\Python36\lib\site-packages\sympy\core\expr.py", line 280, in __float__
raise TypeError("can't convert expression to float")
TypeError: can't convert expression to float
Ich bin kein Informatiker und kann auch sonst wenig und schlecht programmieren. Bitte bitte helft mir. Vielen Dank.
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])
.