SciPy(非线性系统)中ODEINT求解器的数组输入信号

Question

我想传递一个数组输入信号给常微分方程求解器'odeint'而不是定义一个函数。很像线性求解器lsim:

tout, yout, xout = lsim(sys, U=input_sig(time), T=time)

这比定义一个输入函数要方便得多，因为它允许生成各种输入信号并对它们进行操作(窗口或其他)

下面是odeint脚本文档中的一个修改后的示例，其中我将input_sig定义为t

的函数

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

def input_sig(t):
return np.sin(2 * np.pi * 50 * t)
def pend(y, t, b, c):
theta, omega = y
dydt = [omega, -b*omega -c*np.sin(theta) -input_sig(t)]
return dydt

b = 0.25
c = 5.0
y0 = [0, 0.0]
t = np.linspace(0, 10, 101)
sol = odeint(pend, y0, t, args=(b, c))
# what i would like to do is
# sol = odeint(pend, y0, t, input_sig(t), args=(b, c))
# where input_sig(t) is an array
plt.plot(t, sol, label=['theta(t)','omega(t)'])
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

有什么想法吗?我不介意使用除odeint之外的其他解决方案，我不知道是否可能……谢谢你！

Answer 1

正如Lutz Lehman评论中提到的，固定步骤求解器完美地完成了这项工作(我尝试了欧拉正向和龙格库塔):

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

#%% solvers
def rungekutta4(func, y0, t, input_sig, args=()):
n = len(t)
y = np.zeros((n, len(y0)))
y[0] = y0
for i in range(n - 1):
h = t[i+1] - t[i]
k1 = func(y[i], t[i],input_sig[i], *args)
k2 = func(y[i] + k1 * h / 2., t[i] + h / 2.,input_sig[i], *args)
k3 = func(y[i] + k2 * h / 2., t[i] + h / 2.,input_sig[i], *args)
k4 = func(y[i] + k3 * h, t[i] + h,input_sig[i], *args)
y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)
return y
def euler_forward(func, y0, t,input_sig, args=()):
n = len(t)
y = np.zeros((n, len(y0)))
y[0] = y0
for i in range(n - 1):
y[i+1] = y[i] + (t[i+1] - t[i]) * np.array(func(y[i], t[i],input_sig[i], *args))
return y
#%% carac pendulum
g = 9.81        # acceleration of gravity
l = 1           # pendulum rope length
k = 0.8         # air resistance coefficient
m = 1           # mass of the pendulum
b = k/m
c = g/l
d = m *l**2
#%% function
def pend(y, t,input_sig, b, c, d):
theta, omega = y
if isinstance(input_sig,(list,np.ndarray,float)):
dydt = [omega, -b * omega - c * np.sin(theta) - input_sig / d]
else:
dydt = [omega, -b * omega - c * np.sin(theta) - input_sig(t) / d]
return np.array(dydt)

#%% input signal carac
f=20
T = 1          # total duration of the simulation
fs = 20000
dt = 1 / fs      # dt
n_pts= fs * T
t = np.arange(0,n_pts)*dt
def input_sig(t):
return np.sin(2 * np.pi * f * t)
input_sig_array = input_sig(t)
#%% solving
y0 = [0,0]
sol_odeint = odeint(pend, y0, t, args=(input_sig,b, c, d))
sol_e = euler_forward(pend, y0, t,input_sig_array, args=(b,c,d))
sol_r = rungekutta4(pend, y0, t,input_sig_array, args=(b,c,d))
plt.subplot(211)
plt.plot(t, sol_odeint, label=['theta(t)', 'omega(t)'])
plt.plot(t, sol_r,'--',label= ['theta(t)', 'omega(t)'])
plt.plot(t, sol_e,':',label= ['theta(t)', 'omega(t)'])
plt.legend(loc='best')
plt.xlim([0,1/f])
plt.grid()
plt.subplot(212)
plt.plot(t, sol_odeint, label=['theta(t)', 'omega(t)'])
plt.plot(t, sol_r,'--',label= ['theta(t)', 'omega(t)'])
plt.plot(t, sol_e,':',label= ['theta(t)', 'omega(t)'])
plt.legend(loc='best')
plt.xlabel('t, s')
plt.grid()
plt.show()