odeint from scipy.integrate in Python giving wrong result?
Question: I am trying to solve the ivp y'=-y-5 * exp(-t) * sin(5 t), y(0)=1, using the
following code:

    
    
    %pylab inline
    %matplotlib inline
    from scipy.integrate import odeint
    
    def mif(t, y):
        return -y-5*exp(-t)*sin(5*t)
    
    tspan = np.arange(0, 3, 0.000001)
    y0 = 1.0
    y_result = odeint(mif, y0, tspan)
    y_result = y_result[:, 0]  # convert the returned 2D array to a 1D array
    plt.figure()
    plt.plot(tspan, y_result)
    plt.show()
    

However, the plot I get is wrong, it does not match what I obtain, say, with
Matlab or Mathematica. It is actually different from the following alternative
integration:

    
    
    from scipy.integrate import ode
    
    # initialize the 4th order Runge-Kutta solver
    solver = ode(mif).set_integrator('dop853')
    
    # initial value
    y0 = 1.0
    solver.set_initial_value(y0, 0)
    
    values = 1000
    t = np.linspace(0.0001, 3, values)
    y = np.zeros(values)
    
    for ii in range(values):
        y[ii] = solver.integrate(t[ii])[0] #z[0]=u
    

which does yield correct result. What am I doing wrong with the odeint?


Answer: The function arguments change between ode and odeint. For odeint you need

    
    
    def mif(y, t):
    

and for ode

    
    
    def mif(t, y):
    

e.g.

    
    
    %pylab inline
    %matplotlib inline
    from scipy.integrate import odeint
    
    def mif(t,y):
        return y
    
    tspan = np.arange(0, 3, 0.000001)
    y0 = 0.0
    y_result = odeint(mif, y0, tspan)
    plt.figure()
    plt.plot(tspan, y_result)
    plt.show()
    

and

    
    
    from scipy.integrate import ode
    
    def mif(y, t):
        return y
    
    # initialize the 4th order Runge-Kutta solver
    solver = ode(mif).set_integrator('dop853')
    
    # initial value
    y0 = 0.000000
    solver.set_initial_value([y0], 0.0)
    
    values = 1000
    t = np.linspace(0.0000001, 3, values)
    y = np.zeros(values)
    
    for ii in range(values):
        y[ii] = solver.integrate(t[ii]) #z[0]=u
    plt.figure()
    plt.plot(t, y)
    plt.show()
    

