[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