[SciPy-User] Problem with ODE

Issa Karambal issa at aims.ac.za
Wed Dec 19 12:31:48 EST 2012


Hi
Thank to your code I figured out the error I've been getting.
The set_f_params function takes tuple as parameter for the integrators
'dopri5' and 'dopri853' unlike 'vode' integrator.
Thank you very much for your help.

On 19 December 2012 15:39, The Helmbolds <helmrp at yahoo.com> wrote:

> 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
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20121219/8891b1ee/attachment.html>


More information about the SciPy-User mailing list