1

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!

Nuopel
  • 21
  • 5
  • You would need to use a fixed-step solver and also would need to ensure that the implied step size is compatible with the stiffness of the equation. – Lutz Lehmann Jan 05 '23 at 21:11

1 Answers1

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()
Nuopel
  • 21
  • 5
  • You should contemplate the use case for this setup. To make sense for my imagination, the time stepping of the control would be fixed independent of the time stepping of the solver. Thus you would need to up-sample the control to the sampling frequency of the solver, or transform the control data to an interpolation, zero-order hold to remain as step function. – Lutz Lehmann Jan 11 '23 at 11:09