Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_ivp/ivp.py: 11%
167 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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 calling signature is ``fun(t, y)``.
183 Here `t` is a scalar, and there are two options for the ndarray `y`:
184 It can either have shape (n,); then `fun` must return array_like with
185 shape (n,). Alternatively, it can have shape (n, k); then `fun`
186 must return an array_like with shape (n, k), i.e., each column
187 corresponds to a single column in `y`. The choice between the two
188 options is determined by `vectorized` argument (see below). The
189 vectorized implementation allows a faster approximation of the Jacobian
190 by finite differences (required for stiff solvers).
191 t_span : 2-member sequence
192 Interval of integration (t0, tf). The solver starts with t=t0 and
193 integrates until it reaches t=tf. Both t0 and tf must be floats
194 or values interpretable by the float conversion function.
195 y0 : array_like, shape (n,)
196 Initial state. For problems in the complex domain, pass `y0` with a
197 complex data type (even if the initial value is purely real).
198 method : string or `OdeSolver`, optional
199 Integration method to use:
201 * 'RK45' (default): Explicit Runge-Kutta method of order 5(4) [1]_.
202 The error is controlled assuming accuracy of the fourth-order
203 method, but steps are taken using the fifth-order accurate
204 formula (local extrapolation is done). A quartic interpolation
205 polynomial is used for the dense output [2]_. Can be applied in
206 the complex domain.
207 * 'RK23': Explicit Runge-Kutta method of order 3(2) [3]_. The error
208 is controlled assuming accuracy of the second-order method, but
209 steps are taken using the third-order accurate formula (local
210 extrapolation is done). A cubic Hermite polynomial is used for the
211 dense output. Can be applied in the complex domain.
212 * 'DOP853': Explicit Runge-Kutta method of order 8 [13]_.
213 Python implementation of the "DOP853" algorithm originally
214 written in Fortran [14]_. A 7-th order interpolation polynomial
215 accurate to 7-th order is used for the dense output.
216 Can be applied in the complex domain.
217 * 'Radau': Implicit Runge-Kutta method of the Radau IIA family of
218 order 5 [4]_. The error is controlled with a third-order accurate
219 embedded formula. A cubic polynomial which satisfies the
220 collocation conditions is used for the dense output.
221 * 'BDF': Implicit multi-step variable-order (1 to 5) method based
222 on a backward differentiation formula for the derivative
223 approximation [5]_. The implementation follows the one described
224 in [6]_. A quasi-constant step scheme is used and accuracy is
225 enhanced using the NDF modification. Can be applied in the
226 complex domain.
227 * 'LSODA': Adams/BDF method with automatic stiffness detection and
228 switching [7]_, [8]_. This is a wrapper of the Fortran solver
229 from ODEPACK.
231 Explicit Runge-Kutta methods ('RK23', 'RK45', 'DOP853') should be used
232 for non-stiff problems and implicit methods ('Radau', 'BDF') for
233 stiff problems [9]_. Among Runge-Kutta methods, 'DOP853' is recommended
234 for solving with high precision (low values of `rtol` and `atol`).
236 If not sure, first try to run 'RK45'. If it makes unusually many
237 iterations, diverges, or fails, your problem is likely to be stiff and
238 you should use 'Radau' or 'BDF'. 'LSODA' can also be a good universal
239 choice, but it might be somewhat less convenient to work with as it
240 wraps old Fortran code.
242 You can also pass an arbitrary class derived from `OdeSolver` which
243 implements the solver.
244 t_eval : array_like or None, optional
245 Times at which to store the computed solution, must be sorted and lie
246 within `t_span`. If None (default), use points selected by the solver.
247 dense_output : bool, optional
248 Whether to compute a continuous solution. Default is False.
249 events : callable, or list of callables, optional
250 Events to track. If None (default), no events will be tracked.
251 Each event occurs at the zeros of a continuous function of time and
252 state. Each function must have the signature ``event(t, y)`` and return
253 a float. The solver will find an accurate value of `t` at which
254 ``event(t, y(t)) = 0`` using a root-finding algorithm. By default, all
255 zeros will be found. The solver looks for a sign change over each step,
256 so if multiple zero crossings occur within one step, events may be
257 missed. Additionally each `event` function might have the following
258 attributes:
260 terminal: bool, optional
261 Whether to terminate integration if this event occurs.
262 Implicitly False if not assigned.
263 direction: float, optional
264 Direction of a zero crossing. If `direction` is positive,
265 `event` will only trigger when going from negative to positive,
266 and vice versa if `direction` is negative. If 0, then either
267 direction will trigger event. Implicitly 0 if not assigned.
269 You can assign attributes like ``event.terminal = True`` to any
270 function in Python.
271 vectorized : bool, optional
272 Whether `fun` is implemented in a vectorized fashion. Default is False.
273 args : tuple, optional
274 Additional arguments to pass to the user-defined functions. If given,
275 the additional arguments are passed to all user-defined functions.
276 So if, for example, `fun` has the signature ``fun(t, y, a, b, c)``,
277 then `jac` (if given) and any event functions must have the same
278 signature, and `args` must be a tuple of length 3.
279 **options
280 Options passed to a chosen solver. All options available for already
281 implemented solvers are listed below.
282 first_step : float or None, optional
283 Initial step size. Default is `None` which means that the algorithm
284 should choose.
285 max_step : float, optional
286 Maximum allowed step size. Default is np.inf, i.e., the step size is not
287 bounded and determined solely by the solver.
288 rtol, atol : float or array_like, optional
289 Relative and absolute tolerances. The solver keeps the local error
290 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
291 relative accuracy (number of correct digits), while `atol` controls
292 absolute accuracy (number of correct decimal places). To achieve the
293 desired `rtol`, set `atol` to be smaller than the smallest value that
294 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
295 allowable error. If `atol` is larger than ``rtol * abs(y)`` the
296 number of correct digits is not guaranteed. Conversely, to achieve the
297 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
298 than `atol`. If components of y have different scales, it might be
299 beneficial to set different `atol` values for different components by
300 passing array_like with shape (n,) for `atol`. Default values are
301 1e-3 for `rtol` and 1e-6 for `atol`.
302 jac : array_like, sparse_matrix, callable or None, optional
303 Jacobian matrix of the right-hand side of the system with respect
304 to y, required by the 'Radau', 'BDF' and 'LSODA' method. The
305 Jacobian matrix has shape (n, n) and its element (i, j) is equal to
306 ``d f_i / d y_j``. There are three ways to define the Jacobian:
308 * If array_like or sparse_matrix, the Jacobian is assumed to
309 be constant. Not supported by 'LSODA'.
310 * If callable, the Jacobian is assumed to depend on both
311 t and y; it will be called as ``jac(t, y)``, as necessary.
312 For 'Radau' and 'BDF' methods, the return value might be a
313 sparse matrix.
314 * If None (default), the Jacobian will be approximated by
315 finite differences.
317 It is generally recommended to provide the Jacobian rather than
318 relying on a finite-difference approximation.
319 jac_sparsity : array_like, sparse matrix or None, optional
320 Defines a sparsity structure of the Jacobian matrix for a finite-
321 difference approximation. Its shape must be (n, n). This argument
322 is ignored if `jac` is not `None`. If the Jacobian has only few
323 non-zero elements in *each* row, providing the sparsity structure
324 will greatly speed up the computations [10]_. A zero entry means that
325 a corresponding element in the Jacobian is always zero. If None
326 (default), the Jacobian is assumed to be dense.
327 Not supported by 'LSODA', see `lband` and `uband` instead.
328 lband, uband : int or None, optional
329 Parameters defining the bandwidth of the Jacobian for the 'LSODA'
330 method, i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``.
331 Default is None. Setting these requires your jac routine to return the
332 Jacobian in the packed format: the returned array must have ``n``
333 columns and ``uband + lband + 1`` rows in which Jacobian diagonals are
334 written. Specifically ``jac_packed[uband + i - j , j] = jac[i, j]``.
335 The same format is used in `scipy.linalg.solve_banded` (check for an
336 illustration). These parameters can be also used with ``jac=None`` to
337 reduce the number of Jacobian elements estimated by finite differences.
338 min_step : float, optional
339 The minimum allowed step size for 'LSODA' method.
340 By default `min_step` is zero.
342 Returns
343 -------
344 Bunch object with the following fields defined:
345 t : ndarray, shape (n_points,)
346 Time points.
347 y : ndarray, shape (n, n_points)
348 Values of the solution at `t`.
349 sol : `OdeSolution` or None
350 Found solution as `OdeSolution` instance; None if `dense_output` was
351 set to False.
352 t_events : list of ndarray or None
353 Contains for each event type a list of arrays at which an event of
354 that type event was detected. None if `events` was None.
355 y_events : list of ndarray or None
356 For each value of `t_events`, the corresponding value of the solution.
357 None if `events` was None.
358 nfev : int
359 Number of evaluations of the right-hand side.
360 njev : int
361 Number of evaluations of the Jacobian.
362 nlu : int
363 Number of LU decompositions.
364 status : int
365 Reason for algorithm termination:
367 * -1: Integration step failed.
368 * 0: The solver successfully reached the end of `tspan`.
369 * 1: A termination event occurred.
371 message : string
372 Human-readable description of the termination reason.
373 success : bool
374 True if the solver reached the interval end or a termination event
375 occurred (``status >= 0``).
377 References
378 ----------
379 .. [1] J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
380 formulae", Journal of Computational and Applied Mathematics, Vol. 6,
381 No. 1, pp. 19-26, 1980.
382 .. [2] L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
383 of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
384 .. [3] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
385 Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
386 .. [4] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
387 Stiff and Differential-Algebraic Problems", Sec. IV.8.
388 .. [5] `Backward Differentiation Formula
389 <https://en.wikipedia.org/wiki/Backward_differentiation_formula>`_
390 on Wikipedia.
391 .. [6] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
392 COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
393 .. [7] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
394 Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
395 pp. 55-64, 1983.
396 .. [8] L. Petzold, "Automatic selection of methods for solving stiff and
397 nonstiff systems of ordinary differential equations", SIAM Journal
398 on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
399 1983.
400 .. [9] `Stiff equation <https://en.wikipedia.org/wiki/Stiff_equation>`_ on
401 Wikipedia.
402 .. [10] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
403 sparse Jacobian matrices", Journal of the Institute of Mathematics
404 and its Applications, 13, pp. 117-120, 1974.
405 .. [11] `Cauchy-Riemann equations
406 <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
407 Wikipedia.
408 .. [12] `Lotka-Volterra equations
409 <https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations>`_
410 on Wikipedia.
411 .. [13] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
412 Equations I: Nonstiff Problems", Sec. II.
413 .. [14] `Page with original Fortran code of DOP853
414 <http://www.unige.ch/~hairer/software.html>`_.
416 Examples
417 --------
418 Basic exponential decay showing automatically chosen time points.
420 >>> import numpy as np
421 >>> from scipy.integrate import solve_ivp
422 >>> def exponential_decay(t, y): return -0.5 * y
423 >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
424 >>> print(sol.t)
425 [ 0. 0.11487653 1.26364188 3.06061781 4.81611105 6.57445806
426 8.33328988 10. ]
427 >>> print(sol.y)
428 [[2. 1.88836035 1.06327177 0.43319312 0.18017253 0.07483045
429 0.03107158 0.01350781]
430 [4. 3.7767207 2.12654355 0.86638624 0.36034507 0.14966091
431 0.06214316 0.02701561]
432 [8. 7.5534414 4.25308709 1.73277247 0.72069014 0.29932181
433 0.12428631 0.05403123]]
435 Specifying points where the solution is desired.
437 >>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
438 ... t_eval=[0, 1, 2, 4, 10])
439 >>> print(sol.t)
440 [ 0 1 2 4 10]
441 >>> print(sol.y)
442 [[2. 1.21305369 0.73534021 0.27066736 0.01350938]
443 [4. 2.42610739 1.47068043 0.54133472 0.02701876]
444 [8. 4.85221478 2.94136085 1.08266944 0.05403753]]
446 Cannon fired upward with terminal event upon impact. The ``terminal`` and
447 ``direction`` fields of an event are applied by monkey patching a function.
448 Here ``y[0]`` is position and ``y[1]`` is velocity. The projectile starts
449 at position 0 with velocity +10. Note that the integration never reaches
450 t=100 because the event is terminal.
452 >>> def upward_cannon(t, y): return [y[1], -0.5]
453 >>> def hit_ground(t, y): return y[0]
454 >>> hit_ground.terminal = True
455 >>> hit_ground.direction = -1
456 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
457 >>> print(sol.t_events)
458 [array([40.])]
459 >>> print(sol.t)
460 [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
461 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
463 Use `dense_output` and `events` to find position, which is 100, at the apex
464 of the cannonball's trajectory. Apex is not defined as terminal, so both
465 apex and hit_ground are found. There is no information at t=20, so the sol
466 attribute is used to evaluate the solution. The sol attribute is returned
467 by setting ``dense_output=True``. Alternatively, the `y_events` attribute
468 can be used to access the solution at the time of the event.
470 >>> def apex(t, y): return y[1]
471 >>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10],
472 ... events=(hit_ground, apex), dense_output=True)
473 >>> print(sol.t_events)
474 [array([40.]), array([20.])]
475 >>> print(sol.t)
476 [0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
477 1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
478 >>> print(sol.sol(sol.t_events[1][0]))
479 [100. 0.]
480 >>> print(sol.y_events)
481 [array([[-5.68434189e-14, -1.00000000e+01]]), array([[1.00000000e+02, 1.77635684e-15]])]
483 As an example of a system with additional parameters, we'll implement
484 the Lotka-Volterra equations [12]_.
486 >>> def lotkavolterra(t, z, a, b, c, d):
487 ... x, y = z
488 ... return [a*x - b*x*y, -c*y + d*x*y]
489 ...
491 We pass in the parameter values a=1.5, b=1, c=3 and d=1 with the `args`
492 argument.
494 >>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),
495 ... dense_output=True)
497 Compute a dense solution and plot it.
499 >>> t = np.linspace(0, 15, 300)
500 >>> z = sol.sol(t)
501 >>> import matplotlib.pyplot as plt
502 >>> plt.plot(t, z.T)
503 >>> plt.xlabel('t')
504 >>> plt.legend(['x', 'y'], shadow=True)
505 >>> plt.title('Lotka-Volterra System')
506 >>> plt.show()
508 """
509 if method not in METHODS and not (
510 inspect.isclass(method) and issubclass(method, OdeSolver)):
511 raise ValueError("`method` must be one of {} or OdeSolver class."
512 .format(METHODS))
514 t0, tf = map(float, t_span)
516 if args is not None:
517 # Wrap the user's fun (and jac, if given) in lambdas to hide the
518 # additional parameters. Pass in the original fun as a keyword
519 # argument to keep it in the scope of the lambda.
520 try:
521 _ = [*(args)]
522 except TypeError as exp:
523 suggestion_tuple = (
524 "Supplied 'args' cannot be unpacked. Please supply `args`"
525 f" as a tuple (e.g. `args=({args},)`)"
526 )
527 raise TypeError(suggestion_tuple) from exp
529 fun = lambda t, x, fun=fun: fun(t, x, *args)
530 jac = options.get('jac')
531 if callable(jac):
532 options['jac'] = lambda t, x: jac(t, x, *args)
534 if t_eval is not None:
535 t_eval = np.asarray(t_eval)
536 if t_eval.ndim != 1:
537 raise ValueError("`t_eval` must be 1-dimensional.")
539 if np.any(t_eval < min(t0, tf)) or np.any(t_eval > max(t0, tf)):
540 raise ValueError("Values in `t_eval` are not within `t_span`.")
542 d = np.diff(t_eval)
543 if tf > t0 and np.any(d <= 0) or tf < t0 and np.any(d >= 0):
544 raise ValueError("Values in `t_eval` are not properly sorted.")
546 if tf > t0:
547 t_eval_i = 0
548 else:
549 # Make order of t_eval decreasing to use np.searchsorted.
550 t_eval = t_eval[::-1]
551 # This will be an upper bound for slices.
552 t_eval_i = t_eval.shape[0]
554 if method in METHODS:
555 method = METHODS[method]
557 solver = method(fun, t0, y0, tf, vectorized=vectorized, **options)
559 if t_eval is None:
560 ts = [t0]
561 ys = [y0]
562 elif t_eval is not None and dense_output:
563 ts = []
564 ti = [t0]
565 ys = []
566 else:
567 ts = []
568 ys = []
570 interpolants = []
572 events, is_terminal, event_dir = prepare_events(events)
574 if events is not None:
575 if args is not None:
576 # Wrap user functions in lambdas to hide the additional parameters.
577 # The original event function is passed as a keyword argument to the
578 # lambda to keep the original function in scope (i.e., avoid the
579 # late binding closure "gotcha").
580 events = [lambda t, x, event=event: event(t, x, *args)
581 for event in events]
582 g = [event(t0, y0) for event in events]
583 t_events = [[] for _ in range(len(events))]
584 y_events = [[] for _ in range(len(events))]
585 else:
586 t_events = None
587 y_events = None
589 status = None
590 while status is None:
591 message = solver.step()
593 if solver.status == 'finished':
594 status = 0
595 elif solver.status == 'failed':
596 status = -1
597 break
599 t_old = solver.t_old
600 t = solver.t
601 y = solver.y
603 if dense_output:
604 sol = solver.dense_output()
605 interpolants.append(sol)
606 else:
607 sol = None
609 if events is not None:
610 g_new = [event(t, y) for event in events]
611 active_events = find_active_events(g, g_new, event_dir)
612 if active_events.size > 0:
613 if sol is None:
614 sol = solver.dense_output()
616 root_indices, roots, terminate = handle_events(
617 sol, events, active_events, is_terminal, t_old, t)
619 for e, te in zip(root_indices, roots):
620 t_events[e].append(te)
621 y_events[e].append(sol(te))
623 if terminate:
624 status = 1
625 t = roots[-1]
626 y = sol(t)
628 g = g_new
630 if t_eval is None:
631 ts.append(t)
632 ys.append(y)
633 else:
634 # The value in t_eval equal to t will be included.
635 if solver.direction > 0:
636 t_eval_i_new = np.searchsorted(t_eval, t, side='right')
637 t_eval_step = t_eval[t_eval_i:t_eval_i_new]
638 else:
639 t_eval_i_new = np.searchsorted(t_eval, t, side='left')
640 # It has to be done with two slice operations, because
641 # you can't slice to 0th element inclusive using backward
642 # slicing.
643 t_eval_step = t_eval[t_eval_i_new:t_eval_i][::-1]
645 if t_eval_step.size > 0:
646 if sol is None:
647 sol = solver.dense_output()
648 ts.append(t_eval_step)
649 ys.append(sol(t_eval_step))
650 t_eval_i = t_eval_i_new
652 if t_eval is not None and dense_output:
653 ti.append(t)
655 message = MESSAGES.get(status, message)
657 if t_events is not None:
658 t_events = [np.asarray(te) for te in t_events]
659 y_events = [np.asarray(ye) for ye in y_events]
661 if t_eval is None:
662 ts = np.array(ts)
663 ys = np.vstack(ys).T
664 elif ts:
665 ts = np.hstack(ts)
666 ys = np.hstack(ys)
668 if dense_output:
669 if t_eval is None:
670 sol = OdeSolution(ts, interpolants)
671 else:
672 sol = OdeSolution(ti, interpolants)
673 else:
674 sol = None
676 return OdeResult(t=ts, y=ys, sol=sol, t_events=t_events, y_events=y_events,
677 nfev=solver.nfev, njev=solver.njev, nlu=solver.nlu,
678 status=status, message=message, success=status >= 0)