[SciPy-User] Problem with ODE

[The Helmbolds]

On 18 Dec 2012 Issa wrote (heavily abbreviated):

    "I am trying to integrate a simple ode which is succesfull when calling
    'vode' integrator. However,
    with 'dopri5' or 'dop853' as integrator, I have the following error: . . "

I have no idea what the error message means, but the following may be of some help.
Suppose we want to solve the following system of ODEs:
        dx/dt = -Dy
        dy/dt = -Ax
for parameters A = 0.01, D = 0.015 and initial
values x(0) = 1000 and y(0) = 1500. We can
put w = (x, y) and proceed as follows:

    print
    print "ODE-CLASS Results using DOPRI5"
    t0 = 0
    x0 = 1000
    y0 = 1500
    w0 = (x0, y0)
    params = (0.01, 0.015)
    def func_2(t, w, args):    # Note reversal: use (t, w, ...) instead of (w, t, ...).
        A, D = args             # Note no '*'.
        x, y = w
        xdot = -D*y
        ydot = -A*x
        return [xdot, ydot]
    def J_2(t, w, args):       # Jacobian. Note order of (t, w , ...) instead of (w, t, ...).
         A, D = args            # Note no '*'.
         x, y = w
         return [[0, -A], [-D, 0]]
    import scipy
    from scipy.integrate import ode

    # Create an "integrator object" of type "ode"
    #   initialized with the above 'func' and Jacobian.
    intobj = ode(func_2, J_2)

    # Set this 'intobj' to use the 'dopri5' integrator,
    #   using its defaults.
    intobj.set_integrator('dopri5')

    # Set the intobj's initial value and
    #    coresponding intial time.
    intobj.set_initial_value(w0, 0.0)
    # Set the parameters to be used in intobj's
    #   'func' and Jacobian evaluations.
    intobj.set_f_params(params)     # Don't use '*params'.
    intobj.set_jac_params(params)   # Don't use '*params'.

    # Set the final time and time-step to
    #   use for reporting the results.
    tf = 10
    tstep = 1
    t0 = 0.0

    # Display the results, accepting default values
    #   for the remaining options.
    print
    while intobj.successful() and intobj.t < tf:
        intobj.integrate(intobj.t + tstep)
        outstr = '%s, %s' %(intobj.t, intobj.y)
        print outstr
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++        

    print
    print "ODE-CLASS Results using DOP853"
    t0 = 0
    x0 = 1000
    y0 = 1500
    w0 = (x0, y0)
    params = (0.01, 0.015)
    def func_2(t, w, args):    # Note reversal: use (t, w, ...) instead of (w, t, ...).
        A, D = args             # Note no '*'.
        x, y = w
        xdot = -D*y
        ydot = -A*x
        return [xdot, ydot]
    def J_2(t, w, args):       # Jacobian. Note order of (t, w , ...) instead of (w, t, ...).
         A, D = args            # Note no '*'.
         x, y = w
         return [[0, -A], [-D, 0]]
    import scipy
    from scipy.integrate import ode

    # Create an "integrator object" of type "ode"
    #   initialized with the above 'func' and Jacobian.
    intobj = ode(func_2, J_2)

    # Set this 'intobj' to use the 'dop853' integrator, 
    #   using its defaults.
    intobj.set_integrator('dop853')

    # Set the intobj's initial value and
    #    coresponding initial time.
    intobj.set_initial_value(w0, 0.0)
    # Set the parameters to be used in intobj's
    #   'func' and Jacobian evaluations.
    intobj.set_f_params(params)     # Don't use '*params'.
    intobj.set_jac_params(params)   # Don't use '*params'.

    # Set the final time and time-step to
    #   use for reporting the results.
    tf = 10
    tstep = 1
    t0 = 0.0

    # Display the results, accepting default values
    #   for the remaining options.
    print
    while intobj.successful() and intobj.t < tf:
        intobj.integrate(intobj.t + tstep)
        outstr = '%s, %s' %(intobj.t, intobj.y)
        print outstr
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
I get essentially the following output from either of these, or with the theoretical solution.

       1.0, [  977.57443843  1490.1122514 ]
       2.0, [  955.29551487  1480.44802244]
       3.0, [  933.15988742  1471.00586346]
       4.0, [  911.1642357   1461.78435811]
       5.0, [  889.30526033  1452.78212316]
       6.0, [  867.57968241  1443.99780825]
       7.0, [  845.98424307  1435.43009571]
       8.0, [  824.51570296  1427.07770039]
       9.0, [  803.17084174  1418.93936939]
       10.0, [  781.94645766  1411.01388197]

Hope this helps.
Bob H