pythonnumpyscipydifferential-equationsnumerical-integration

How to use dorpi5 or dop853 in Python


I have looked through scipy.integrate.ode but I am unable to find out how to actually use these integration methods, dorpi5 and dop853.

I would like to try integrating ode integration python versus mathematica my python code with these two methods to see how it effects the results but don't know how.


Solution

  • You call the method set_integrator on the ode class with either 'dopri5' or 'dop853' as its argument.

    Here's an example:

    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.integrate import ode
    
    
    def fun(t, z, omega):
        """
        Right hand side of the differential equations
          dx/dt = -omega * y
          dy/dt = omega * x
        """
        x, y = z
        f = [-omega*y, omega*x]
        return f
    
    # Create an `ode` instance to solve the system of differential
    # equations defined by `fun`, and set the solver method to 'dop853'.
    solver = ode(fun)
    solver.set_integrator('dop853')
    
    # Give the value of omega to the solver. This is passed to
    # `fun` when the solver calls it.
    omega = 2 * np.pi
    solver.set_f_params(omega)
    
    # Set the initial value z(0) = z0.
    t0 = 0.0
    z0 = [1, -0.25]
    solver.set_initial_value(z0, t0)
    
    # Create the array `t` of time values at which to compute
    # the solution, and create an array to hold the solution.
    # Put the initial value in the solution array.
    t1 = 2.5
    N = 75
    t = np.linspace(t0, t1, N)
    sol = np.empty((N, 2))
    sol[0] = z0
    
    # Repeatedly call the `integrate` method to advance the
    # solution to time t[k], and save the solution in sol[k].
    k = 1
    while solver.successful() and solver.t < t1:
        solver.integrate(t[k])
        sol[k] = solver.y
        k += 1
    
    # Plot the solution...
    plt.plot(t, sol[:,0], label='x')
    plt.plot(t, sol[:,1], label='y')
    plt.xlabel('t')
    plt.grid(True)
    plt.legend()
    plt.show()
    

    This generates the following plot:

    plot