0
$\begingroup$

My time step with the function scipy.integrate.solve_ivp is not decreasing in t_span fluctuating (reaching values below or higher than t_end) which is a problem for my application. Can anyone tell me if, first, it is a normal situation and, second, if it is possible to make the time step continuously decrease in t_span?

Here is an example of the time step that the function is passing through

-32.0
-24.572060540292647 <--- This first drop is the problem
-31.895622059576212
-31.843433089364318
-31.791244119152424
-31.395622059576212
-31.197811029788106
-31.0
t_span = [-32, -31]

Here the part of my code using this function. I use THETA to compute X and depending on the value of X, I know if I am in the compression stroke, combustion phase, or expansion stroke but this time drop put me in the combustion phase but I am actually in the compression stroke.

def rates( THETA, Y):
   YPRIME = zeros(NY);
   VOL, X , EM = auxiliary(THETA);
   print(THETA)
   M = EM*MNOT;
   VOL, DV, Ac = volume(THETA)
   AA = (DV + VOL*BLOWBY/OMEGA)/M;
   C1 = HEAT*(Ac)/OMEGA/M/1000;
   C0 = sqrt(X);
   P = Y[0];
   TB = Y[1];
   TU = Y[2];

   if ( X > 0.999 ):
       #  EXPANSION
       ierr, YB, HL, xxx, xxx, VB, xxx, CP, xxx, DVDT, DVDP = ecp( TB, P, PHI, fuel_type );
       if ( ierr != 0 ):
           print('Error in ECP({0}, {1}, {2}): {3}\n'.format(TB, P, PHI, ierr))                      
    
       BB = C1/CP*DVDT*TB*(1-TW/TB);
       CC = 0;
       DD = 1/CP*TB*DVDT**2 + DVDP;
       EE = 0;
    
       YPRIME[0] = (AA + BB + CC)/(DD + EE);
       YPRIME[1] = -C1/CP*(TB-TW) + 1/CP*DVDT*TB*YPRIME[0];
       YPRIME[2] = 0;
    
   elif ( X > 0.001 ):
       #  COMBUSTION
       xxx, HU, xxx, xxx, VU, xxx, CPU, xxx, DVDTU, DVDPU = farg( TU, P, PHI, F, fuel_type );
       ierr, YB, HB, xxx, xxx, VB, xxx, CPB, xxx, DVDTB, DVDPB = ecp( TB, P, PHI, fuel_type );
       if ( ierr != 0 ):
           print('Error in ECP({0}, {1}, {2}): {3}\n'.format(TB, P, PHI, ierr))      
       BB = C1*(1/CPB*TB*DVDTB*C0*(1-TW/TB) + 1/CPU*TU*DVDTU*(1-C0)*(1-TW/TU));
       DX = xb(THETA,THETAS,THETAB) [1]
       CC = -(VB-VU)*DX - DVDTB*(HU-HB)/CPB*(DX-(X-X**2)*BLOWBY/OMEGA);
       DD = X*(VB**2/CPB/TB*(TB/VB*DVDTB)**2 + DVDPB);
       EE = (1-X)*(1/CPU*TU*DVDTU**2 + DVDPU);
       HL = (1-X**2)*HU + X**2*HB;
    
       YPRIME[0] = (AA + BB + CC)/(DD + EE);
       YPRIME[1] = -C1/CPB/C0*(TB-TW) + 1/CPB*TB*DVDTB*YPRIME[0] + (HU-HB)/CPB*(DX/X - (1-X)*BLOWBY/OMEGA);
       YPRIME[2] = -C1/CPU/(1+C0)*(TU-TW) + 1/CPU*TU*DVDTU*YPRIME[0];
    
   else:
       #  COMPRESSION            
       xxx, HL, xxx, xxx, xxx, xxx, CP, xxx, DVDT, DVDP = farg( TU, P, PHI, F, fuel_type )
    
       BB = C1*1/CP*TU*DVDT*(1-TW/TU)
       CC = 0
       DD = 0
       EE = 1/CP*TU*DVDT**2 + DVDP
       YPRIME[0] = ( AA + BB + CC )/(DD + EE)
       YPRIME[1] = 0
       YPRIME[2] = -C1/CP*(TU-TW) + 1/CP*TU*DVDT*YPRIME[0]

   # 1/omega is s/rad, so convert to s/deg
   for JJ in range(NY) : 
       YPRIME[JJ] = YPRIME[JJ]*pi/180

   return YPRIME

def integrate( THETA, THETAE, Y, X):
   sol = scipy.integrate.solve_ivp(rates, array([THETA, THETAE]), Y, method='RK23', args=(X,))
   return sol
$\endgroup$

1 Answer 1

1
$\begingroup$

Scipy's solve_ivp has indeed an issue where it may sometimes go beyond the final time at the first step. Since you do not provide an initial step size (via the first_step argument), it attempts computing one based on a simple method (described e.g. in Hairer ang Wanner, "Solving Ordinary Differential Equations I", sect II.4). See the source code here: https://github.com/scipy/scipy/blob/f43f724827c9e2d40a04b7ab4bd2d46a2632c241/scipy/integrate/_ivp/common.py#L68

This routine is however such that the proposed initial step size may be larger than the time interval to be covered by the integration. While waiting for this to be fixed, definitely set a first_step value, e.g. (tmax-tmin)/10 or something like that.

Note that the LSODA wrapper from solve_ivp does not suffer from this issue, since it relies on a Fortran code that has its own initial time step determination method, which ensures that it remains below tmax-tmin.

For more info, see e.g. https://github.com/scipy/scipy/pull/17348 .

$\endgroup$
3
  • $\begingroup$ I really wish there were good ODE solver packages in Python that aren't just wrappers to other implemented methods. $\endgroup$
    – whpowell96
    Apr 10 at 14:39
  • 1
    $\begingroup$ The solvers available in Scipy are in my opinion quite good for prototyping and small enough systems. Wrapping other methods seems rather a.sensible thing to do for me. There are a few missing bits of functionality here and there but overall I find its current state sufficient. $\endgroup$
    – Laurent90
    Apr 10 at 15:43
  • $\begingroup$ Yeah LSODA and SUNDIALS wrappers can handle most things but it's nice to have more methods and optimization options within the same interface as in Julia DiffEq.jl or PETSC. $\endgroup$
    – whpowell96
    Apr 10 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.