6

I wanted to store the different integration steps taken by the solver itself when I call it :

solver1.integrate(t_end)

So I did a while loop and enabled the step option setting its value to True:

while solver1.successful() and solver1.t < t0+dt:
    solver1.integrate(t_end,step=True)
    time.append(solver1.t)

Then I plot y, the result of integration and here comes my issue. I have instabilities which appear in a located area :

y enabling the step option of the solver

I thought it was because of the loop or something like that so I checked the result removing the step :

while solver1.successful() and solver1.t < t0+dt:
    solver1.integrate(t_end)

And surprise ... I have the correct result :

y disabling the step option of the solver

It's a quite weird situation ... I'd be grateful if someone of your guys could help me with this issue.

EDIT :

To set the solver I do :

solver1 = ode(y_dot,jac).set_integrator('vode',with_jacobian=True)
solver1.set_initial_value(x0,t0)

And I store the result using .append()

kuider
  • 103
  • 5
  • Can you show some more of your code, how you set up the solver and store the result for plotting? – silvado Mar 20 '13 at 15:36
  • Of course, I've just edited my question. – kuider Mar 21 '13 at 08:28
  • You still haven't shown how you are actually storing the current ODE state which you are plotting, presuming that the plots show one of the ODE state variables. – Nikolas Apr 02 '13 at 18:41

1 Answers1

2

When you set step=True you are indirectly giving the vode._integrator.runner (a Fortran subroutine) an instruction to use itask=2, and the default is itask=1. You can get more details about this runner doing:

r._integrator.runner?

In SciPy 0.12.0 documentation you will not find an explanation about what is going on for the different itask=1 or itask=2, but you can find it here:

ITASK  = An index specifying the task to be performed.
!          Input only. ITASK has the following values and meanings.
!          1  means normal computation of output values of y(t) at
!             t = TOUT(by overshooting and interpolating).
!          2  means take one step only and return.
!          3  means stop at the first internal mesh point at or
!             beyond t = TOUT and return.
!          4  means normal computation of output values of y(t) at
!             t = TOUT but without overshooting t = TCRIT.
!             TCRIT must be input as RUSER(1). TCRIT may be equal to
!             or beyond TOUT, but not behind it in the direction of
!             integration. This option is useful if the problem
!             has a singularity at or beyond t = TCRIT.
!          5  means take one step, without passing TCRIT, and return.
!             TCRIT must be input as RUSER(1).
Saullo G. P. Castro
  • 56,802
  • 26
  • 179
  • 234