Python3 ode solver with non-zero initial conditions failing

Question

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?

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

Answer 1

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