How to use dorpi5 or dop853 in Python

dustin picture dustin · Apr 26, 2013 · Viewed 7.9k times · Source

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.

Answer

Warren Weckesser picture Warren Weckesser · Apr 26, 2013

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