How to correctly implement a time dependent variable when using scipy.integrate.odeint

Question

I'm trying to solve a system of odes basically resembling this one but with one more spring and damper ==> http://scipy-cookbook.readthedocs.io/items/CoupledSpringMassSystem.html

I have a slight problem, though, because one of the parameters I want to implement is time dependent. My first attempt is the following one :

import scipy as sci
import numpy as np
import matplotlib.pyplot as plt

def bump(t):    
    if t <= (0.25 / 6.9):
        return (0.075 * (1 - np.cos(np.pi * 8 * 6.9 * t)))
    else:
        return 0 


def membre_droite(w, t, p):
    x1,y1,x2,y2,x3,y3 = w
    m1,m2,m3,k1,k2,k3,l1,l2,l3,c2,c3 = p
    f = [y1,
         (-k1 * (x1 - l1 - bump(t)) + k2 * (x2 - x1 - l2) + c2 * (y2 - y1)) / m1,
          y2,
          (-c2 * (y2 - y1) - k2 * (x2 - x1 - l2) + k3 * (x3 - x2 - l3) + c3 * (y3 - y2)) / m2,
          y3,
          (-c3 * (y3 - y2) - k3 * (x3 - x2 - l3)) / m3]
    return f



# Initial values
x11 = 0.08
y11 = 0
x22 = 0.35
y22 = 0
x33 = 0.6
y33 = 0

# Parameters
m1 = 90
m2 = 4000
m3 = 105
k1 = 250000
k2 = 25000
k3 = 30000
l1 = 0.08
l2 = x22-x11
l3 = x33-x22
c2 = 2500
c3 = 850

# Initial paramaters regrouped + time array 
time=np.linspace(0.0, 5, 1001)
w0 = [x11,y11,x22,y22,x33,y33]
p0 = [m1,m2,m3,k1,k2,k3,l1,l2,l3,c2,c3]
x1,y1,x2,y2,x3,y3 = sci.integrate.odeint(membre_droite, w0, time, args=(p0,)).T

plt.plot(time,x1,'b')
plt.plot(time,x2,'g')
plt.plot(time,x3,'r') 
plt.plot(time,y2,'yellow')
plt.plot(time,y3,'black')                                      
plt.xlabel('t')
plt.grid(True)
plt.legend((r'$x_1$', r'$x_2$', r'$x_3$', r'$y_2$', r'$y_3$'))
plt.show()

The error I get is :

if t <= (0.25 / 6.9):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

I've looked for similar cases and I came across this topic ==> Solving a system of odes (with changing constant!) using scipy.integrate.odeint?

I've then attempted to adapt my code to this format :

import scipy as sci
import numpy as np
import matplotlib.pyplot as plt

def bump(t):    
    if t <= (0.25 / 6.9):
        return (0.075 * (1 - np.cos(np.pi * 8 * 6.9 * t)))
    else:
        return 0


def membre_droite(w, t, bump):
    x1,y1,x2,y2,x3,y3 = w
    f = [y1,
         (-250000 * (x1 - x11 - bump(t)) + 25000 * (x2 - x1 - x22 + x11) + 2500 * (y2-y1)) / 90,
         y2,
         (-2500 * (y2 - y1) - 25000 * (x2 - x1 - x22 + x11) + 30000 * (x3 - x2 - x33 + x22) + 850 * (y3 - y2)) / 4000,
         y3,
         (-850 * (y3 - y2) - 30000 * (x3 - x2 - x33 + x22)) / 105]
    return f



# Initial values
x11 = 0.08
y11 = 0
x22 = 0.35
y22 = 0
x33 = 0.6
y33 = 0



# Initial paramaters regrouped + time array 
time = np.linspace(0.0, 5, 1001)
w0 = [x11,y11,x22,y22,x33,y33]

x1,y1,x2,y2,x3,y3 = sci.integrate.odeint(membre_droite, w0, time, args=(bump,)).T

plt.plot(time,x1,'b')
plt.plot(time,x2,'g')
plt.plot(time,x3,'r') 
plt.plot(time,y2,'yellow')
plt.plot(time,y3,'black')                                      
plt.xlabel('t')
plt.grid(True)
plt.legend((r'$x_1$', r'$x_2$', r'$x_3$', r'$y_2$', r'$y_3$'))
plt.show()

Reading the previous link, it should have worked but I get another error :

(-250000 * (x1 - x11 - bump(t)) + 25000 * (x2 - x1 - x22 + x11) + 2500 * (y2 - y1)) / 90,
TypeError: 'list' object is not callable

I have tested your code with some modifications like: `from scipy.integrate import odeint` and `x1,y1,x2,y2,x3,y3 = odeint(membre_droite, w0, time, args=(bump,)).T` and it works fine, what version of python and libraries are you using. — eyllanesc, May 24 '17 at 16:57
I'm using the 3.6 version. What do you mean for the libraries, though? (I'm relatively new to all this so I don't see what I could have changed about that) — Evogaiser, May 24 '17 at 17:02
Well actually for some reason when I switched to "from scipy.integrate import odeint" like you said, it worked. I don't know what could have caused that but thank you very much. — Evogaiser, May 24 '17 at 17:10
I'm going to put that answer and mark it as the answer please. — eyllanesc, May 24 '17 at 17:11

score 4 · Accepted Answer · answered May 24 '17 at 17:12

Change to:

from scipy.integrate import odeint

and

x1,y1,x2,y2,x3,y3 = odeint(membre_droite, w0, time, args=(bump,)).T

Complete code:

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

def bump(t):    
    if t <= (0.25 / 6.9):
        return (0.075 * (1 - np.cos(np.pi * 8 * 6.9 * t)))
    else:
        return 0

def membre_droite(w, t, bump):
    x1,y1,x2,y2,x3,y3 = w
    f = [y1,
         (-250000 * (x1 - x11 - bump(t)) + 25000 * (x2 - x1 - x22 + x11) + 2500 * (y2 - y1)) / 90,
         y2,
         (-2500 * (y2 - y1) - 25000 * (x2 - x1 - x22 + x11) + 30000 * (x3 - x2 - x33 + x22) + 850 * (y3 - y2)) / 4000,
         y3,
         (-850 * (y3 - y2) - 30000 * (x3 - x2 - x33 + x22)) / 105]
    return f



# Initial values
x11 = 0.08
y11 = 0
x22 = 0.35
y22 = 0
x33 = 0.6
y33 = 0



# Initial paramaters regrouped + time array 
time = np.linspace(0.0, 5, 1001)
w0 = [x11,y11,x22,y22,x33,y33]

x1,y1,x2,y2,x3,y3 = odeint(membre_droite, w0, time, args=(bump,)).T

plt.plot(time,x1,'b')
plt.plot(time,x2,'g')
plt.plot(time,x3,'r') 
plt.plot(time,y2,'yellow')
plt.plot(time,y3,'black')                                      
plt.xlabel('t')
plt.grid(True)
plt.legend((r'$x_1$', r'$x_2$', r'$x_3$', r'$y_2$', r'$y_3$'))
plt.show()

Max Herrmann · Answer 2 · 2023-06-30T19:49:32.547

Introduce a time-dependent variable u by setting its value at each time step, e.g.

u(t_0) = 1, u(t_1) = 1, u(t_2) = 0, u(t_3) = 0, ...

Then interpolate between these values using numpy.interp().

Here is a minimal example:

from scipy.integrate import odeint
import numpy as np

# time vector
N = 128
t = np.linspace(0, 4, N)

# time-dependent variable u = [1, 1, ..., 1, 0, ..., 0]
u = np.array([0]*N)
u[np.nonzero(t < 0.5)] = 1

# right-hand side of ode
def rhs(x, t, u, tt):

    # find input value
    u_val = np.interp(t, tt, u, left=0, right=0)
    
    # dx = -x + u
    return -x + u_val

# initial value
x0 = 0

# solve
sol = odeint(rhs, x0, t, args=(u, t))

How to correctly implement a time dependent variable when using scipy.integrate.odeint

2 Answers2