Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_ivp/ivp.py: 11%
168 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import inspect
2import numpy as np
3from .bdf import BDF
4from .radau import Radau
5from .rk import RK23, RK45, DOP853
6from .lsoda import LSODA
7from scipy.optimize import OptimizeResult
8from .common import EPS, OdeSolution
9from .base import OdeSolver
12METHODS = {'RK23': RK23,
13 'RK45': RK45,
14 'DOP853': DOP853,
15 'Radau': Radau,
16 'BDF': BDF,
17 'LSODA': LSODA}
20MESSAGES = {0: "The solver successfully reached the end of the integration interval.",
21 1: "A termination event occurred."}
24class OdeResult(OptimizeResult):
25 pass
28def prepare_events(events):
29 """Standardize event functions and extract is_terminal and direction."""
30 if callable(events):
31 events = (events,)
33 if events is not None:
34 is_terminal = np.empty(len(events), dtype=bool)
35 direction = np.empty(len(events))
36 for i, event in enumerate(events):
37 try:
38 is_terminal[i] = event.terminal
39 except AttributeError:
40 is_terminal[i] = False
42 try:
43 direction[i] = event.direction
44 except AttributeError:
45 direction[i] = 0
46 else:
47 is_terminal = None
48 direction = None
50 return events, is_terminal, direction
53def solve_event_equation(event, sol, t_old, t):
54 """Solve an equation corresponding to an ODE event.
56 The equation is ``event(t, y(t)) = 0``, here ``y(t)`` is known from an
57 ODE solver using some sort of interpolation. It is solved by
58 `scipy.optimize.brentq` with xtol=atol=4*EPS.
60 Parameters
61 ----------
62 event : callable
63 Function ``event(t, y)``.
64 sol : callable
65 Function ``sol(t)`` which evaluates an ODE solution between `t_old`
66 and `t`.
67 t_old, t : float
68 Previous and new values of time. They will be used as a bracketing
69 interval.
71 Returns
72 -------
73 root : float
74 Found solution.
75 """
76 from scipy.optimize import brentq
77 return brentq(lambda t: event(t, sol(t)), t_old, t,
78 xtol=4 * EPS, rtol=4 * EPS)
81def handle_events(sol, events, active_events, is_terminal, t_old, t):
82 """Helper function to handle events.
84 Parameters
85 ----------
86 sol : DenseOutput
87 Function ``sol(t)`` which evaluates an ODE solution between `t_old`
88 and `t`.
89 events : list of callables, length n_events
90 Event functions with signatures ``event(t, y)``.
91 active_events : ndarray
92 Indices of events which occurred.
93 is_terminal : ndarray, shape (n_events,)
94 Which events are terminal.
95 t_old, t : float
96 Previous and new values of time.
98 Returns
99 -------
100 root_indices : ndarray
101 Indices of events which take zero between `t_old` and `t` and before
102 a possible termination.
103 roots : ndarray
104 Values of t at which events occurred.
105 terminate : bool
106 Whether a terminal event occurred.
107 """
108 roots = [solve_event_equation(events[event_index], sol, t_old, t)
109 for event_index in active_events]
111 roots = np.asarray(roots)
113 if np.any(is_terminal[active_events]):
114 if t > t_old:
115 order = np.argsort(roots)
116 else:
117 order = np.argsort(-roots)
118 active_events = active_events[order]
119 roots = roots[order]
120 t = np.nonzero(is_terminal[active_events])[0][0]
121 active_events = active_events[:t + 1]
122 roots = roots[:t + 1]
123 terminate = True
124 else:
125 terminate = False
127 return active_events, roots, terminate
130def find_active_events(g, g_new, direction):
131 """Find which event occurred during an integration step.
133 Parameters
134 ----------
135 g, g_new : array_like, shape (n_events,)
136 Values of event functions at a current and next points.
137 direction : ndarray, shape (n_events,)
138 Event "direction" according to the definition in `solve_ivp`.
140 Returns
141 -------
142 active_events : ndarray
143 Indices of events which occurred during the step.
144 """
145 g, g_new = np.asarray(g), np.asarray(g_new)
146 up = (g <= 0) & (g_new >= 0)
147 down = (g >= 0) & (g_new <= 0)
148 either = up | down
149 mask = (up & (direction > 0) |
150 down & (direction < 0) |
151 either & (direction == 0))
153 return np.nonzero(mask)[0]
156def solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False,
157 events=None, vectorized=False, args=None, **options):
158 """Solve an initial value problem for a system of ODEs.
160 This function numerically integrates a system of ordinary differential
161 equations given an initial value::
163 dy / dt = f(t, y)
164 y(t0) = y0
166 Here t is a 1-D independent variable (time), y(t) is an
167 N-D vector-valued function (state), and an N-D
168 vector-valued function f(t, y) determines the differential equations.
169 The goal is to find y(t) approximately satisfying the differential
170 equations, given an initial value y(t0)=y0.
172 Some of the solvers support integration in the complex domain, but note
173 that for stiff ODE solvers, the right-hand side must be
174 complex-differentiable (satisfy Cauchy-Riemann equations [11]_).
175 To solve a problem in the complex domain, pass y0 with a complex data type.
176 Another option always available is to rewrite your problem for real and
177 imaginary parts separately.
179 Parameters
180 ----------
181 fun : callable
182 Right-hand side of the system: the time derivative of the state ``y``
183 at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a
184 scalar and ``y`` is an ndarray with ``len(y) = len(y0)``. ``fun`` must
185 return an array of the same shape as ``y``. See `vectorized` for more
186 information.
187 t_span : 2-member sequence
188 Interval of integration (t0, tf). The solver starts with t=t0 and
189 integrates until it reaches t=tf. Both t0 and tf must be floats
190 or values interpretable by the float conversion function.
191 y0 : array_like, shape (n,)
192 Initial state. For problems in the complex domain, pass `y0` with a
193 complex data type (even if the initial value is purely real).
194 method : string or `OdeSolver`, optional
195 Integration method to use:
197 * 'RK45' (default): Explicit Runge-Kutta method of order 5(4) [1]_.
198 The error is controlled assuming accuracy of the fourth-order
199 method, but steps are taken using the fifth-order accurate
200 formula (local extrapolation is done). A quartic interpolation
201 polynomial is used for the dense output [2]_. Can be applied in
202 the complex domain.
203 * 'RK23': Explicit Runge-Kutta method of order 3(2) [3]_. The error
204 is controlled assuming accuracy of the second-order method, but
205 steps are taken using the third-order accurate formula (local
206 extrapolation is done). A cubic Hermite polynomial is used for the
207 dense output. Can be applied in the complex domain.
208 * 'DOP853': Explicit Runge-Kutta method of order 8 [13]_.
209 Python implementation of the "DOP853" algorithm originally
210 written in Fortran [14]_. A 7-th order interpolation polynomial
211 accurate to 7-th order is used for the dense output.
212 Can be applied in the complex domain.
213 * 'Radau': Implicit Runge-Kutta method of the Radau IIA family of
214 order 5 [4]_. The error is controlled with a third-order accurate
215 embedded formula. A cubic polynomial which satisfies the
216 collocation conditions is used for the dense output.
217 * 'BDF': Implicit multi-step variable-order (1 to 5) method based
218 on a backward differentiation formula for the derivative
219 approximation [5]_. The implementation follows the one described
220 in [6]_. A quasi-constant step scheme is used and accuracy is
221 enhanced using the NDF modification. Can be applied in the
222 complex domain.
223 * 'LSODA': Adams/BDF method with automatic stiffness detection and
224 switching [7]_, [8]_. This is a wrapper of the Fortran solver
225 from ODEPACK.
227 Explicit Runge-Kutta methods ('RK23', 'RK45', 'DOP853') should be used
228 for non-stiff problems and implicit methods ('Radau', 'BDF') for
229 stiff problems [9]_. Among Runge-Kutta methods, 'DOP853' is recommended
230 for solving with high precision (low values of `rtol` and `atol`).
232 If not sure, first try to run 'RK45'. If it makes unusually many
233 iterations, diverges, or fails, your problem is likely to be stiff and
234 you should use 'Radau' or 'BDF'. 'LSODA' can also be a good universal
235 choice, but it might be somewhat less convenient to work with as it
236 wraps old Fortran code.
238 You can also pass an arbitrary class derived from `OdeSolver` which
239 implements the solver.
240 t_eval : array_like or None, optional
241 Times at which to store the computed solution, must be sorted and lie
242 within `t_span`. If None (default), use points selected by the solver.
243 dense_output : bool, optional
244 Whether to compute a continuous solution. Default is False.
245 events : callable, or list of callables, optional
246 Events to track. If None (default), no events will be tracked.
247 Each event occurs at the zeros of a continuous function of time and
248 state. Each function must have the signature ``event(t, y)`` and return
249 a float. The solver will find an accurate value of `t` at which
250 ``event(t, y(t)) = 0`` using a root-finding algorithm. By default, all
251 zeros will be found. The solver looks for a sign change over each step,
252 so if multiple zero crossings occur within one step, events may be
253 missed. Additionally each `event` function might have the following
254 attributes:
256 terminal: bool, optional
257 Whether to terminate integration if this event occurs.
258 Implicitly False if not assigned.
259 direction: float, optional
260 Direction of a zero crossing. If `direction` is positive,
261 `event` will only trigger when going from negative to positive,
262 and vice versa if `direction` is negative. If 0, then either
263 direction will trigger event. Implicitly 0 if not assigned.
265 You can assign attributes like ``event.terminal = True`` to any
266 function in Python.
267 vectorized : bool, optional
268 Whether `fun` can be called in a vectorized fashion. Default is False.
270 If ``vectorized`` is False, `fun` will always be called with ``y`` of
271 shape ``(n,)``, where ``n = len(y0)``.
273 If ``vectorized`` is True, `fun` may be called with ``y`` of shape
274 ``(n, k)``, where ``k`` is an integer. In this case, `fun` must behave
275 such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of
276 the returned array is the time derivative of the state corresponding
277 with a column of ``y``).
279 Setting ``vectorized=True`` allows for faster finite difference
280 approximation of the Jacobian by methods 'Radau' and 'BDF', but
281 will result in slower execution for other methods and for 'Radau' and
282 'BDF' in some circumstances (e.g. small ``len(y0)``).
283 args : tuple, optional
284 Additional arguments to pass to the user-defined functions. If given,
285 the additional arguments are passed to all user-defined functions.
286 So if, for example, `fun` has the signature ``fun(t, y, a, b, c)``,
287 then `jac` (if given) and any event functions must have the same
288 signature, and `args` must be a tuple of length 3.
289 **options
290 Options passed to a chosen solver. All options available for already
291 implemented solvers are listed below.
292 first_step : float or None, optional
293 Initial step size. Default is `None` which means that the algorithm
294 should choose.
295 max_step : float, optional
296 Maximum allowed step size. Default is np.inf, i.e., the step size is not
297 bounded and determined solely by the solver.
298 rtol, atol : float or array_like, optional
299 Relative and absolute tolerances. The solver keeps the local error
300 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
301 relative accuracy (number of correct digits), while `atol` controls
302 absolute accuracy (number of correct decimal places). To achieve the
303 desired `rtol`, set `atol` to be smaller than the smallest value that
304 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
305 allowable error. If `atol` is larger than ``rtol * abs(y)`` the
306 number of correct digits is not guaranteed. Conversely, to achieve the
307 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
308 than `atol`. If components of y have different scales, it might be
309 beneficial to set different `atol` values for different components by
310 passing array_like with shape (n,) for `atol`. Default values are
311 1e-3 for `rtol` and 1e-6 for `atol`.
312 jac : array_like, sparse_matrix, callable or None, optional
313 Jacobian matrix of the right-hand side of the system with respect
314 to y, required by the 'Radau', 'BDF' and 'LSODA' method. The
315 Jacobian matrix has shape (n, n) and its element (i, j) is equal to
316 ``d f_i / d y_j``. There are three ways to define the Jacobian:
318 * If array_like or sparse_matrix, the Jacobian is assumed to
319 be constant. Not supported by 'LSODA'.
320 * If callable, the Jacobian is assumed to depend on both
321 t and y; it will be called as ``jac(t, y)``, as necessary.
322 For 'Radau' and 'BDF' methods, the return value might be a
323 sparse matrix.
324 * If None (default), the Jacobian will be approximated by
325 finite differences.
327 It is generally recommended to provide the Jacobian rather than
328 relying on a finite-difference approximation.
329 jac_sparsity : array_like, sparse matrix or None, optional
330 Defines a sparsity structure of the Jacobian matrix for a finite-
331 difference approximation. Its shape must be (n, n). This argument
332 is ignored if `jac` is not `None`. If the Jacobian has only few
333 non-zero elements in *each* row, providing the sparsity structure
334 will greatly speed up the computations [10]_. A zero entry means that
335 a corresponding element in the Jacobian is always zero. If None
336 (default), the Jacobian is assumed to be dense.
337 Not supported by 'LSODA', see `lband` and `uband` instead.
338 lband, uband : int or None, optional
339 Parameters defining the bandwidth of the Jacobian for the 'LSODA'
340 method, i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``.
341 Default is None. Setting these requires your jac routine to return the
342 Jacobian in the packed format: the returned array must have ``n``
343 columns and ``uband + lband + 1`` rows in which Jacobian diagonals are
344 written. Specifically ``jac_packed[uband + i - j , j] = jac[i, j]``.
345 The same format is used in `scipy.linalg.solve_banded` (check for an
346 illustration). These parameters can be also used with ``jac=None`` to
347 reduce the number of Jacobian elements estimated by finite differences.
348 min_step : float, optional
349 The minimum allowed step size for 'LSODA' method.
350 By default `min_step` is zero.
352 Returns
353 -------
354 Bunch object with the following fields defined:
355 t : ndarray, shape (n_points,)
356 Time points.
357 y : ndarray, shape (n, n_points)
358 Values of the solution at `t`.
359 sol : `OdeSolution` or None
360 Found solution as `OdeSolution` instance; None if `dense_output` was
361 set to False.
362 t_events : list of ndarray or None
363 Contains for each event type a list of arrays at which an event of
364 that type event was detected. None if `events` was None.
365 y_events : list of ndarray or None
366 For each value of `t_events`, the corresponding value of the solution.
367 None if `events` was None.
368 nfev : int
369 Number of evaluations of the right-hand side.
370 njev : int
371 Number of evaluations of the Jacobian.
372 nlu : int
373 Number of LU decompositions.
374 status : int
375 Reason for algorithm termination:
377 * -1: Integration step failed.
378 * 0: The solver successfully reached the end of `tspan`.
379 * 1: A termination event occurred.
381 message : string
382 Human-readable description of the termination reason.
383 success : bool
384 True if the solver reached the interval end or a termination event
385 occurred (``status >= 0``).
387 References
388 ----------
389 .. [1] J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
390 formulae", Journal of Computational and Applied Mathematics, Vol. 6,
391 No. 1, pp. 19-26, 1980.
392 .. [2] L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
393 of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
394 .. [3] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
395 Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
396 .. [4] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
397 Stiff and Differential-Algebraic Problems", Sec. IV.8.
398 .. [5] `Backward Differentiation Formula
399 <https://en.wikipedia.org/wiki/Backward_differentiation_formula>`_
400 on Wikipedia.
401 .. [6] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
402 COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
403 .. [7] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
404 Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
405 pp. 55-64, 1983.
406 .. [8] L. Petzold, "Automatic selection of methods for solving stiff and
407 nonstiff systems of ordinary differential equations", SIAM Journal
408 on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
409 1983.
410 .. [9] `Stiff equation <https://en.wikipedia.org/wiki/Stiff_equation>`_ on
411 Wikipedia.
412 .. [10] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
413 sparse Jacobian matrices", Journal of the Institute of Mathematics
414 and its Applications, 13, pp. 117-120, 1974.
415 .. [11] `Cauchy-Riemann equations
416 <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
417 Wikipedia.
418 .. [12] `Lotka-Volterra equations
419 <https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations>`_
420 on Wikipedia.
421 .. [13] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
422 Equations I: Nonstiff Problems", Sec. II.
423 .. [14] `Page with original Fortran code of DOP853
424 <http://www.unige.ch/~hairer/software.html>`_.
426 Examples
427 --------
428 Basic exponential decay showing automatically chosen time points.
430 >>> import numpy as np
431 >>> from scipy.integrate import solve_ivp
432 >>> def exponential_decay(t, y): return -0.5 * y
433 >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
434 >>> print(sol.t)
435 [ 0. 0.11487653 1.26364188 3.06061781 4.81611105 6.57445806
436 8.33328988 10. ]
437 >>> print(sol.y)
438 [[2. 1.88836035 1.06327177 0.43319312 0.18017253 0.07483045
439 0.03107158 0.01350781]
440 [4. 3.7767207 2.12654355 0.86638624 0.36034507 0.14966091
441 0.06214316 0.02701561]
442 [8. 7.5534414 4.25308709 1.73277247 0.72069014 0.29932181
443 0.12428631 0.05403123]]
445 Specifying points where the solution is desired.
447 >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
448 ... t_eval=[0, 1, 2, 4, 10])
449 >>> print(sol.t)
450 [ 0 1 2 4 10]
451 >>> print(sol.y)
452 [[2. 1.21305369 0.73534021 0.27066736 0.01350938]
453 [4. 2.42610739 1.47068043 0.54133472 0.02701876]
454 [8. 4.85221478 2.94136085 1.08266944 0.05403753]]
456 Cannon fired upward with terminal event upon impact. The ``terminal`` and
457 ``direction`` fields of an event are applied by monkey patching a function.
458 Here ``y[0]`` is position and ``y[1]`` is velocity. The projectile starts
459 at position 0 with velocity +10. Note that the integration never reaches
460 t=100 because the event is terminal.
462 >>> def upward_cannon(t, y): return [y[1], -0.5]
463 >>> def hit_ground(t, y): return y[0]
464 >>> hit_ground.terminal = True
465 >>> hit_ground.direction = -1
466 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
467 >>> print(sol.t_events)
468 [array([40.])]
469 >>> print(sol.t)
470 [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
471 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
473 Use `dense_output` and `events` to find position, which is 100, at the apex
474 of the cannonball's trajectory. Apex is not defined as terminal, so both
475 apex and hit_ground are found. There is no information at t=20, so the sol
476 attribute is used to evaluate the solution. The sol attribute is returned
477 by setting ``dense_output=True``. Alternatively, the `y_events` attribute
478 can be used to access the solution at the time of the event.
480 >>> def apex(t, y): return y[1]
481 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10],
482 ... events=(hit_ground, apex), dense_output=True)
483 >>> print(sol.t_events)
484 [array([40.]), array([20.])]
485 >>> print(sol.t)
486 [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
487 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
488 >>> print(sol.sol(sol.t_events[1][0]))
489 [100. 0.]
490 >>> print(sol.y_events)
491 [array([[-5.68434189e-14, -1.00000000e+01]]), array([[1.00000000e+02, 1.77635684e-15]])]
493 As an example of a system with additional parameters, we'll implement
494 the Lotka-Volterra equations [12]_.
496 >>> def lotkavolterra(t, z, a, b, c, d):
497 ... x, y = z
498 ... return [a*x - b*x*y, -c*y + d*x*y]
499 ...
501 We pass in the parameter values a=1.5, b=1, c=3 and d=1 with the `args`
502 argument.
504 >>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),
505 ... dense_output=True)
507 Compute a dense solution and plot it.
509 >>> t = np.linspace(0, 15, 300)
510 >>> z = sol.sol(t)
511 >>> import matplotlib.pyplot as plt
512 >>> plt.plot(t, z.T)
513 >>> plt.xlabel('t')
514 >>> plt.legend(['x', 'y'], shadow=True)
515 >>> plt.title('Lotka-Volterra System')
516 >>> plt.show()
518 """
519 if method not in METHODS and not (
520 inspect.isclass(method) and issubclass(method, OdeSolver)):
521 raise ValueError("`method` must be one of {} or OdeSolver class."
522 .format(METHODS))
524 t0, tf = map(float, t_span)
526 if args is not None:
527 # Wrap the user's fun (and jac, if given) in lambdas to hide the
528 # additional parameters. Pass in the original fun as a keyword
529 # argument to keep it in the scope of the lambda.
530 try:
531 _ = [*(args)]
532 except TypeError as exp:
533 suggestion_tuple = (
534 "Supplied 'args' cannot be unpacked. Please supply `args`"
535 f" as a tuple (e.g. `args=({args},)`)"
536 )
537 raise TypeError(suggestion_tuple) from exp
539 def fun(t, x, fun=fun):
540 return fun(t, x, *args)
541 jac = options.get('jac')
542 if callable(jac):
543 options['jac'] = lambda t, x: jac(t, x, *args)
545 if t_eval is not None:
546 t_eval = np.asarray(t_eval)
547 if t_eval.ndim != 1:
548 raise ValueError("`t_eval` must be 1-dimensional.")
550 if np.any(t_eval < min(t0, tf)) or np.any(t_eval > max(t0, tf)):
551 raise ValueError("Values in `t_eval` are not within `t_span`.")
553 d = np.diff(t_eval)
554 if tf > t0 and np.any(d <= 0) or tf < t0 and np.any(d >= 0):
555 raise ValueError("Values in `t_eval` are not properly sorted.")
557 if tf > t0:
558 t_eval_i = 0
559 else:
560 # Make order of t_eval decreasing to use np.searchsorted.
561 t_eval = t_eval[::-1]
562 # This will be an upper bound for slices.
563 t_eval_i = t_eval.shape[0]
565 if method in METHODS:
566 method = METHODS[method]
568 solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
570 if t_eval is None:
571 ts = [t0]
572 ys = [y0]
573 elif t_eval is not None and dense_output:
574 ts = []
575 ti = [t0]
576 ys = []
577 else:
578 ts = []
579 ys = []
581 interpolants = []
583 events, is_terminal, event_dir = prepare_events(events)
585 if events is not None:
586 if args is not None:
587 # Wrap user functions in lambdas to hide the additional parameters.
588 # The original event function is passed as a keyword argument to the
589 # lambda to keep the original function in scope (i.e., avoid the
590 # late binding closure "gotcha").
591 events = [lambda t, x, event=event: event(t, x, *args)
592 for event in events]
593 g = [event(t0, y0) for event in events]
594 t_events = [[] for _ in range(len(events))]
595 y_events = [[] for _ in range(len(events))]
596 else:
597 t_events = None
598 y_events = None
600 status = None
601 while status is None:
602 message = solver.step()
604 if solver.status == 'finished':
605 status = 0
606 elif solver.status == 'failed':
607 status = -1
608 break
610 t_old = solver.t_old
611 t = solver.t
612 y = solver.y
614 if dense_output:
615 sol = solver.dense_output()
616 interpolants.append(sol)
617 else:
618 sol = None
620 if events is not None:
621 g_new = [event(t, y) for event in events]
622 active_events = find_active_events(g, g_new, event_dir)
623 if active_events.size > 0:
624 if sol is None:
625 sol = solver.dense_output()
627 root_indices, roots, terminate = handle_events(
628 sol, events, active_events, is_terminal, t_old, t)
630 for e, te in zip(root_indices, roots):
631 t_events[e].append(te)
632 y_events[e].append(sol(te))
634 if terminate:
635 status = 1
636 t = roots[-1]
637 y = sol(t)
639 g = g_new
641 if t_eval is None:
642 ts.append(t)
643 ys.append(y)
644 else:
645 # The value in t_eval equal to t will be included.
646 if solver.direction > 0:
647 t_eval_i_new = np.searchsorted(t_eval, t, side='right')
648 t_eval_step = t_eval[t_eval_i:t_eval_i_new]
649 else:
650 t_eval_i_new = np.searchsorted(t_eval, t, side='left')
651 # It has to be done with two slice operations, because
652 # you can't slice to 0th element inclusive using backward
653 # slicing.
654 t_eval_step = t_eval[t_eval_i_new:t_eval_i][::-1]
656 if t_eval_step.size > 0:
657 if sol is None:
658 sol = solver.dense_output()
659 ts.append(t_eval_step)
660 ys.append(sol(t_eval_step))
661 t_eval_i = t_eval_i_new
663 if t_eval is not None and dense_output:
664 ti.append(t)
666 message = MESSAGES.get(status, message)
668 if t_events is not None:
669 t_events = [np.asarray(te) for te in t_events]
670 y_events = [np.asarray(ye) for ye in y_events]
672 if t_eval is None:
673 ts = np.array(ts)
674 ys = np.vstack(ys).T
675 elif ts:
676 ts = np.hstack(ts)
677 ys = np.hstack(ys)
679 if dense_output:
680 if t_eval is None:
681 sol = OdeSolution(ts, interpolants)
682 else:
683 sol = OdeSolution(ti, interpolants)
684 else:
685 sol = None
687 return OdeResult(t=ts, y=ys, sol=sol, t_events=t_events, y_events=y_events,
688 nfev=solver.nfev, njev=solver.njev, nlu=solver.nlu,
689 status=status, message=message, success=status >= 0)