scipy odeint: Array einfügen

[blueme]

Hallo zusammen,

mein Problem ist das Folgende:
ich würde gerne eine DGL 2ter Ordnung lösen, allerdings hat diese DGL ein Störglied, welches mir als eine Liste vorliegt.
Diese DGL würde ich gerne mittel scipy.odeint lösen, allerdings schaffe ich nicht die Liste des Störgliedes in die Gleichung einzubauen. Wenn ich eine einzelne Zahl oder np.cos(t) nehme funktioniert es, nur eben bei einer Liste nicht.

Code: Alles auswählen

def model(z,t):
    u = np.array([1,1,1,1,1,1,1,1,1,1])
    #u = np.cos(t)
    return np.array([z[1], u -6*z[1]-10*z[0]])

u ist eine Beispiel Liste die ich gerne zu "return" hinzufügen würde, folgende Fehlermeldung kommt dann allerdings: "TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'"

Ich hoffe ich konnte das Problem verständlich beschreiben.
Über Hilfe wäre ich sehr dankbar

Hier unten noch der dazugehörige Beispielcode.

Code: Alles auswählen

import matplotlib.pyplot as plt
from scipy.integrate import odeint
import numpy as np

# function that returns dz/dt
def model(z,t):
    u = np.array([1,1,1,1,1,1,1,1,1,1])
    #u = np.cos(t)
    return np.array([z[1], u -6*z[1]-10*z[0]])

# initial condition
z0 = [0,4]

# time points
t = np.linspace(0,5,10)

# solve ODE
z = odeint(model,z0,t)

# plot results
plt.plot(t,z[:,0],'b-',label= "x")
plt.ylabel('response')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

[einfachTobi]

Zunächst solltest du `solve_ivp`statt `odeint` verwenden (siehe Kasten hier).
Man beachte die geänderte Reihenfolge in der Modellfunktion:

Code: Alles auswählen

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import numpy as np

def model(t, z):
    #u = np.array([1,1,1,1,1,1,1,1,1,1])
    u = np.cos(t)
    return [z[1], u-6*z[1]-10*z[0]]

z0 = [0, 4]
result = solve_ivp(model, (0, 5), z0, max_step=0.1)

print(result)

plt.plot(result.t, result.y[0], label='result[0]')
plt.plot(result.t, result.y[1], label='result[1]')
plt.ylabel('response')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

Wie kann ich mir das Störglied in Wirklichkeit vorstellen anstatt einer Liste mit lauter 1en? Wie ist deine Modellvorstellung davon? Wann soll welches Element der Liste des Störgliedes verwendet werden? Kannst du das Originalproblem näher erläutern?

[blueme]

Mein Originalproblem ist folgendes: ich habe eine Eingangsgröße (das u), die zu bestimmten Zeitpunkten abgefragt wurde, diese Zeitpunkte habe ich auch als Liste. Desweiteren habe ich die Ausgangsgröße, welche zu den selben Zeitpunkten abgefragt wird. Mein Ziel ist es die DGL zu lösen und die Vorfaktoren ( in unserem Beispiel noch 6 und 10, aber eigentlich variabel) so anzupassen, dass die Lösung der DGL gleich der Ausgangsgröße ist (ungefähr).
Aktuell habe ich das Problem gelöst indem in die Liste interpoliert habe (linear) und so konnte ich sie einfügen.
Was ist der Vorteil von `solve_ivp`statt `odeint` ?

Vielen Dank für deine Hilfe!

[einfachTobi]

Vorteile: z. B. Kompatibilität mit kommenden scipy-Versionen, Auswahlmöglichkeit aus mehreren Solvern, Verwendung von Abbruchkriterien, Rückmeldung über Erfolg/Abbruch.
Das klingt für mich nach Curve-Fitting. Lässt sich also leicht automatisieren:

Code: Alles auswählen

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.optimize import curve_fit
import numpy as np

def model(t, z):
    u = np.cos(t)
    return [z[1], u-6*z[1]-10*z[0]]

def model_fit(t, z, a, b):
    u = np.cos(t)
    return [z[1], u-a*z[1]-b*z[0]]

def fitfunc(t, a, b):
    def my_ode(t, z):
        u = np.cos(t)
        return [z[1], u-a*z[1]-b*z[0]]
    ode_solution = solve_ivp(my_ode, (0, 1), z0, t_eval=tspan)
    return ode_solution.y[0]

tspan = [0, 0.2, 0.4, 0.6, 0.8, 1]
data = [0.0081,  0.5,  0.53,  0.42,  0.32,  0.23]
z0 = [0, 4]

# konventionell lösen mit den Beispielwerten 6 und 10
result = solve_ivp(model, (0, 1), z0, max_step=0.01)

# Curve-Fitting, um die Parameter a und b von my_ode an die Messwerte anzupassen
k_fit, kcov = curve_fit(fitfunc, tspan, data)
print("k_fit:", k_fit)

# Lösung der DGL mit den ermittelten Werten a und b zur Kontrolle
fit = solve_ivp(model_fit, t_span=(0, 1), y0=z0, args=(k_fit[0], k_fit[1]), max_step=0.01)

plt.plot(tspan, data, 'ro', label='Messwerte')
plt.plot(result.t, result.y[0], label='6 und 10')
plt.plot(fit.t, fit.y[0], label='fit mit a und b')
plt.ylabel('response')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()