Array input signal for ODEINT solver in SciPy (nonlinear system)

Question

I want to pass an array input signal to ordinary differential equations solver 'odeint' instead of defining a function. Much like the linear solver lsim :

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

Which is much more convenient than defining a input function as it allows to generate various input signal and operate on them (windowing or whatever)

In the following is a modified example from the odeint scipy documentation where i define input_sig as a function of 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()

Any ideas ? i don't mind using other solvers than odeint, i don't know if it's possible... Thank You!

Answer 1

As mentioned in Lutz Lehman comment a fixed-step solver did the job perfectly (I tried Euler forward and Runge Kutta):

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()