How to correctly implement forced mass spring system ode solver with a force time dependent variable using scipy.integrate.odeint

Question

I am trying, for each of the 5 cases, to integrate numerically through the odeint function, a spring mass function, the parameter F being variable with time. However, it presents the error:

line 244, in odeint int(bool(tfirst)))

ValueError: setting an array element with a sequence.

from scipy.integrate import odeint as ode
import matplotlib.pyplot as graph
import numpy as np

def spring(y,t,par):

    m = par[0]
    k = par[1]
    c = par[2]
    F = par[3] 

    ydot=[0,0]
    ydot[0] = y[1]
    F = np.array(F)
    ydot[1] = F/m-c*y[1]/m-k*y[0]/m
    return ydot

#cases:[ti,tf,x0,xdot0,c]
A=[0.0,0.3,0.1,0.0,37.7]
B=[0.0,3.0,0.0,1.0,37.7]
C=[0.0,3.0,0.1,1.0,377]
D=[0.0,5.0,0.0,0.0,37.7]
E=[0.0,3.0,0.0,0.0,37.7]

cases = [A, B, C, D, E] #0,1,2,3,4

m=10.0
k=3553.0
par = [m,k]
h = 0.01
cont = 0

for case in cases: 
    x0 = case[2]
    xdot0 = case[3]
    y = [x0,xdot0]
    par.append(case[4])
    ti = case[0]
    tf = case[1]
    t = np.arange(ti, tf,h)

    F = []
    for time in t:
        if cont == 3:
            F.append(1000*np.sin(np.pi*time+np.pi/2))
        elif (cont == 4) and (time >= 0.5): 
            F.append(1000)
        else:
            F.append(0)
    cont = cont + 1
    par.append(F) #F is a vector

    Y = ode(spring,y,t,args=(par,))

    Xn = Y[:,0]
    Vn = Y[:,1]

    graph.plot(t,Xn)
    graph.show()

    graph.plot(t,Vn)
    graph.show()

    graph.plot(Xn,Vn)
    graph.show()

You do not clean the parameter list between cases, how do you think the solver will extract the current force at time `t` from the list passed as parameter? Especially as that is not the task of the solver but of the derivatives function. `spring` does not do anything in that regard. — Lutz Lehmann, Mar 25 '20 at 08:16

score 0 · Answer 1 · answered Mar 25 '20 at 00:05

There are other problems in the code which will arise when you try to run, after having fixed that problem. But just in the interest of being precise, I'll answer just your question.

The problem is in this line of code:

ydot[1] = F/m-c*y[1]/m-k*y[0]/m

The F here is currently fed in as a list. In maths if you divide a vector by a scalar it is assumed that you perform element-wise division. Python lists don't replicate this behaviour, but numpy arrays do. So just convert your list to numpy array like so:

F = np.array(F)
ydot[1] = F/m-c*y[1]/m-k*y[0]/m

DrBwts · Answer 2 · 2020-03-25T14:29:53.447

0

The problem is with the return value from spring(y,t,par)

odeint is expecting 2 single values for ydot[0] & ydot[1] to be returned from spring() but is instead getting a single value in ydot[0] and then an array of length len(F) in ydot[1].

Change,

ydot[1] = F / m - c * y[1] / m - k * y[0] / m

to

ydot[1] = F[0] / m - c * y[1] / m - k * y[0] / m

and check what happens..

edited Mar 25 '20 at 14:29

answered Mar 25 '20 at 14:21

DrBwts

3,470
6
38
62

score 0 · Answer 3 · answered Mar 25 '20 at 15:07

You are completely wrong in passing the value array F to the derivatives function without any means of relating its entries to the time t that the function is evaluated at. Remember that the arguments y and t are the state vector and (single, scalar) time that the solver needs to compute the next stage in the numerical method. The easiest way is to pass F as function so that its value is can be directly computed.

def spring(y,t,par):

    m,k,c,F = par 
    return [ y[1], F(t)/m-c*y[1]/m-k*y[0]/m ]

#cases:[ti,tf,x0,xdot0,c]
A=[0.0,0.3,0.1,0.0,37.7]
B=[0.0,3.0,0.0,1.0,37.7]
C=[0.0,3.0,0.1,1.0,377]
D=[0.0,5.0,0.0,0.0,37.7]
E=[0.0,3.0,0.0,0.0,37.7]

cases = [A, B, C, D, E] #0,1,2,3,4

m=10.0
k=3553.0
h = 0.01

It is more pythonic to use tuple assignments to unpack parameter tuples, doing in one line what was previously distributed over several, also making the order of parameters better visible.

Remove some unnecessary complications, such as having an extra counting mechanism, use the enumerate(list) mechanism. Also, do not use an explicitly constructed parameter list par when it is much easier to construct it ad-hoc in the args optional argument (this would be different if there were hundreds of systematically constructed entries).

In a general setting F might be an interpolation function generated with methods from scipy.interpolate. Here it is sufficient to transform the given equations into lambda expressions to define the corresponding scalar functions.

for cont,case in enumerate(cases): 
    ti,tf,x0,xdot0,c = case 
    y = [x0,xdot0]
    t = np.arange(ti, tf, h)

    F = lambda t: 0
    if cont == 3:
            F = lambda t: 1000*np.sin(np.pi*t+np.pi/2)
    elif (cont == 4): 
            F = lambda t:  1000 if t >= 0.5 else 0

    Y = ode(spring,y,t,args=([m,k,c,F],))

    Xn,Vn = Y.T

    graph.figure(figsize=(9,3))
    graph.subplot(1,3,1); graph.plot(t,Xn)
    graph.subplot(1,3,2); graph.plot(t,Vn)

    graph.subplot(1,3,3); graph.plot(Xn,Vn)
    graph.tight_layout(); graph.show()

Everything else proceeds as before, using sub-plots to assemble the different graphs in one picture. The 4th example for instance gives the chaotically oscillating result

How to correctly implement forced mass spring system ode solver with a force time dependent variable using scipy.integrate.odeint

3 Answers3