0

I am trying to simulate a mass spring system with a single degree of freedom. For the time integration I am using the ode function from scipy. Also, I am comparing the numerics with the analytical solution.

Running my script results in the same frequency response as the analytical solution, but there is a jump in the amplitude after the first time step. I cannot find the source of error in my code. Does anyone have an idea?

enter image description here

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


def main():

    # spring stiffness, mass, initial displacement
    K, M, u0 = 1.

    v0 = 0.

    # max time and time step size
    T = 10.
    dt = 1E-5

    u, v, t = ODE(K,M,T,dt,u0,v0)

    # analytical solution
    u_analytical = u0 * np.cos(np.sqrt(K/M)*t)

    plt.plot(t,u_analytical)
    plt.plot(t,u)
    plt.legend(['analytic', 'numeric'])
    plt.show()



def MSD(t,states,K,M):

    u, v = np.reshape(states, (2, -1))

    a = -K/M*u

    states[0] = v
    states[1] = a

    return states


def ODE(K,M,T,dt,u0,v0):

    states0 = np.array([u0, v0])

    r = ode(MSD)

    r.set_integrator('vode', method='bdf', order=5, nsteps=1000, rtol=1E-8)

    r.set_initial_value(states0, 0)
    r.set_f_params(K,M)

    time_array = np.linspace(0, T, T/dt + 1)
    u = [u0]
    v = [v0]

    for t in range(1, len(time_array)):

        states = r.integrate(time_array[t])

        u.append(states[0])
        v.append(states[1])

        assert r.successful()

    u = np.asarray(u)
    v = np.asarray(v)

    return u, v, time_array

main()

2 Answers 2

0

Your problem is that you try to re-use the state vector. Now if the integrator does something like in this simple Euler step,

y += f(t,y)*dt

then changing y as side effect of the call of f will dramatically change the result.


If you assume that you have possibly multiple instances of the system as indicated by the first line

    u, v = np.reshape(states, (2, -1))

then you should care about flattening the resulting derivatives vector

    return np.concatenate([v,a])
2
  • Thanks for the explanation. As stated below I found out that using a list instead of a numpy array also solves the problem. Still, I don't really understand the issue. Why can't I just fill the array 'states' with indexing as in my example? P.S. I accepted your answer
    – The Dens
    May 4, 2020 at 18:54
  • Because python passes arguments as references. With you re-using state you might and obviously do overwrite an internal variable of the solver. A similar trick that tries to avoid the memory allocation of the derivatives vector is to make a class that has the return vector as member and the ODE function as method that returns always the same vector. This might slightly speed up larger computations (but for those better use code-writing methods like jitcode). May 4, 2020 at 22:09
0

I found the solution to my problem.

When the function MBD returns a list instead of a numpy array it works:

def MSD_fixed(t,states,K,M):

    u, v = np.reshape(states, (2, -1))

    a = -K/M*u

    return [v, a]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.