Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_ivp/rk.py: 32%
191 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
2from .base import OdeSolver, DenseOutput
3from .common import (validate_max_step, validate_tol, select_initial_step,
4 norm, warn_extraneous, validate_first_step)
5from . import dop853_coefficients
7# Multiply steps computed from asymptotic behaviour of errors by this.
8SAFETY = 0.9
10MIN_FACTOR = 0.2 # Minimum allowed decrease in a step size.
11MAX_FACTOR = 10 # Maximum allowed increase in a step size.
14def rk_step(fun, t, y, f, h, A, B, C, K):
15 """Perform a single Runge-Kutta step.
17 This function computes a prediction of an explicit Runge-Kutta method and
18 also estimates the error of a less accurate method.
20 Notation for Butcher tableau is as in [1]_.
22 Parameters
23 ----------
24 fun : callable
25 Right-hand side of the system.
26 t : float
27 Current time.
28 y : ndarray, shape (n,)
29 Current state.
30 f : ndarray, shape (n,)
31 Current value of the derivative, i.e., ``fun(x, y)``.
32 h : float
33 Step to use.
34 A : ndarray, shape (n_stages, n_stages)
35 Coefficients for combining previous RK stages to compute the next
36 stage. For explicit methods the coefficients at and above the main
37 diagonal are zeros.
38 B : ndarray, shape (n_stages,)
39 Coefficients for combining RK stages for computing the final
40 prediction.
41 C : ndarray, shape (n_stages,)
42 Coefficients for incrementing time for consecutive RK stages.
43 The value for the first stage is always zero.
44 K : ndarray, shape (n_stages + 1, n)
45 Storage array for putting RK stages here. Stages are stored in rows.
46 The last row is a linear combination of the previous rows with
47 coefficients
49 Returns
50 -------
51 y_new : ndarray, shape (n,)
52 Solution at t + h computed with a higher accuracy.
53 f_new : ndarray, shape (n,)
54 Derivative ``fun(t + h, y_new)``.
56 References
57 ----------
58 .. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
59 Equations I: Nonstiff Problems", Sec. II.4.
60 """
61 K[0] = f
62 for s, (a, c) in enumerate(zip(A[1:], C[1:]), start=1):
63 dy = np.dot(K[:s].T, a[:s]) * h
64 K[s] = fun(t + c * h, y + dy)
66 y_new = y + h * np.dot(K[:-1].T, B)
67 f_new = fun(t + h, y_new)
69 K[-1] = f_new
71 return y_new, f_new
74class RungeKutta(OdeSolver):
75 """Base class for explicit Runge-Kutta methods."""
76 C: np.ndarray = NotImplemented
77 A: np.ndarray = NotImplemented
78 B: np.ndarray = NotImplemented
79 E: np.ndarray = NotImplemented
80 P: np.ndarray = NotImplemented
81 order: int = NotImplemented
82 error_estimator_order: int = NotImplemented
83 n_stages: int = NotImplemented
85 def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,
86 rtol=1e-3, atol=1e-6, vectorized=False,
87 first_step=None, **extraneous):
88 warn_extraneous(extraneous)
89 super().__init__(fun, t0, y0, t_bound, vectorized,
90 support_complex=True)
91 self.y_old = None
92 self.max_step = validate_max_step(max_step)
93 self.rtol, self.atol = validate_tol(rtol, atol, self.n)
94 self.f = self.fun(self.t, self.y)
95 if first_step is None:
96 self.h_abs = select_initial_step(
97 self.fun, self.t, self.y, self.f, self.direction,
98 self.error_estimator_order, self.rtol, self.atol)
99 else:
100 self.h_abs = validate_first_step(first_step, t0, t_bound)
101 self.K = np.empty((self.n_stages + 1, self.n), dtype=self.y.dtype)
102 self.error_exponent = -1 / (self.error_estimator_order + 1)
103 self.h_previous = None
105 def _estimate_error(self, K, h):
106 return np.dot(K.T, self.E) * h
108 def _estimate_error_norm(self, K, h, scale):
109 return norm(self._estimate_error(K, h) / scale)
111 def _step_impl(self):
112 t = self.t
113 y = self.y
115 max_step = self.max_step
116 rtol = self.rtol
117 atol = self.atol
119 min_step = 10 * np.abs(np.nextafter(t, self.direction * np.inf) - t)
121 if self.h_abs > max_step:
122 h_abs = max_step
123 elif self.h_abs < min_step:
124 h_abs = min_step
125 else:
126 h_abs = self.h_abs
128 step_accepted = False
129 step_rejected = False
131 while not step_accepted:
132 if h_abs < min_step:
133 return False, self.TOO_SMALL_STEP
135 h = h_abs * self.direction
136 t_new = t + h
138 if self.direction * (t_new - self.t_bound) > 0:
139 t_new = self.t_bound
141 h = t_new - t
142 h_abs = np.abs(h)
144 y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
145 self.B, self.C, self.K)
146 scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol
147 error_norm = self._estimate_error_norm(self.K, h, scale)
149 if error_norm < 1:
150 if error_norm == 0:
151 factor = MAX_FACTOR
152 else:
153 factor = min(MAX_FACTOR,
154 SAFETY * error_norm ** self.error_exponent)
156 if step_rejected:
157 factor = min(1, factor)
159 h_abs *= factor
161 step_accepted = True
162 else:
163 h_abs *= max(MIN_FACTOR,
164 SAFETY * error_norm ** self.error_exponent)
165 step_rejected = True
167 self.h_previous = h
168 self.y_old = y
170 self.t = t_new
171 self.y = y_new
173 self.h_abs = h_abs
174 self.f = f_new
176 return True, None
178 def _dense_output_impl(self):
179 Q = self.K.T.dot(self.P)
180 return RkDenseOutput(self.t_old, self.t, self.y_old, Q)
183class RK23(RungeKutta):
184 """Explicit Runge-Kutta method of order 3(2).
186 This uses the Bogacki-Shampine pair of formulas [1]_. The error is controlled
187 assuming accuracy of the second-order method, but steps are taken using the
188 third-order accurate formula (local extrapolation is done). A cubic Hermite
189 polynomial is used for the dense output.
191 Can be applied in the complex domain.
193 Parameters
194 ----------
195 fun : callable
196 Right-hand side of the system: the time derivative of the state ``y``
197 at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a
198 scalar and ``y`` is an ndarray with ``len(y) = len(y0)``. ``fun`` must
199 return an array of the same shape as ``y``. See `vectorized` for more
200 information.
201 t0 : float
202 Initial time.
203 y0 : array_like, shape (n,)
204 Initial state.
205 t_bound : float
206 Boundary time - the integration won't continue beyond it. It also
207 determines the direction of the integration.
208 first_step : float or None, optional
209 Initial step size. Default is ``None`` which means that the algorithm
210 should choose.
211 max_step : float, optional
212 Maximum allowed step size. Default is np.inf, i.e., the step size is not
213 bounded and determined solely by the solver.
214 rtol, atol : float and array_like, optional
215 Relative and absolute tolerances. The solver keeps the local error
216 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
217 relative accuracy (number of correct digits), while `atol` controls
218 absolute accuracy (number of correct decimal places). To achieve the
219 desired `rtol`, set `atol` to be smaller than the smallest value that
220 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
221 allowable error. If `atol` is larger than ``rtol * abs(y)`` the
222 number of correct digits is not guaranteed. Conversely, to achieve the
223 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
224 than `atol`. If components of y have different scales, it might be
225 beneficial to set different `atol` values for different components by
226 passing array_like with shape (n,) for `atol`. Default values are
227 1e-3 for `rtol` and 1e-6 for `atol`.
228 vectorized : bool, optional
229 Whether `fun` may be called in a vectorized fashion. False (default)
230 is recommended for this solver.
232 If ``vectorized`` is False, `fun` will always be called with ``y`` of
233 shape ``(n,)``, where ``n = len(y0)``.
235 If ``vectorized`` is True, `fun` may be called with ``y`` of shape
236 ``(n, k)``, where ``k`` is an integer. In this case, `fun` must behave
237 such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of
238 the returned array is the time derivative of the state corresponding
239 with a column of ``y``).
241 Setting ``vectorized=True`` allows for faster finite difference
242 approximation of the Jacobian by methods 'Radau' and 'BDF', but
243 will result in slower execution for this solver.
245 Attributes
246 ----------
247 n : int
248 Number of equations.
249 status : string
250 Current status of the solver: 'running', 'finished' or 'failed'.
251 t_bound : float
252 Boundary time.
253 direction : float
254 Integration direction: +1 or -1.
255 t : float
256 Current time.
257 y : ndarray
258 Current state.
259 t_old : float
260 Previous time. None if no steps were made yet.
261 step_size : float
262 Size of the last successful step. None if no steps were made yet.
263 nfev : int
264 Number evaluations of the system's right-hand side.
265 njev : int
266 Number of evaluations of the Jacobian. Is always 0 for this solver as it does not use the Jacobian.
267 nlu : int
268 Number of LU decompositions. Is always 0 for this solver.
270 References
271 ----------
272 .. [1] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
273 Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
274 """
275 order = 3
276 error_estimator_order = 2
277 n_stages = 3
278 C = np.array([0, 1/2, 3/4])
279 A = np.array([
280 [0, 0, 0],
281 [1/2, 0, 0],
282 [0, 3/4, 0]
283 ])
284 B = np.array([2/9, 1/3, 4/9])
285 E = np.array([5/72, -1/12, -1/9, 1/8])
286 P = np.array([[1, -4 / 3, 5 / 9],
287 [0, 1, -2/3],
288 [0, 4/3, -8/9],
289 [0, -1, 1]])
292class RK45(RungeKutta):
293 """Explicit Runge-Kutta method of order 5(4).
295 This uses the Dormand-Prince pair of formulas [1]_. The error is controlled
296 assuming accuracy of the fourth-order method accuracy, but steps are taken
297 using the fifth-order accurate formula (local extrapolation is done).
298 A quartic interpolation polynomial is used for the dense output [2]_.
300 Can be applied in the complex domain.
302 Parameters
303 ----------
304 fun : callable
305 Right-hand side of the system. The calling signature is ``fun(t, y)``.
306 Here ``t`` is a scalar, and there are two options for the ndarray ``y``:
307 It can either have shape (n,); then ``fun`` must return array_like with
308 shape (n,). Alternatively it can have shape (n, k); then ``fun``
309 must return an array_like with shape (n, k), i.e., each column
310 corresponds to a single column in ``y``. The choice between the two
311 options is determined by `vectorized` argument (see below).
312 t0 : float
313 Initial time.
314 y0 : array_like, shape (n,)
315 Initial state.
316 t_bound : float
317 Boundary time - the integration won't continue beyond it. It also
318 determines the direction of the integration.
319 first_step : float or None, optional
320 Initial step size. Default is ``None`` which means that the algorithm
321 should choose.
322 max_step : float, optional
323 Maximum allowed step size. Default is np.inf, i.e., the step size is not
324 bounded and determined solely by the solver.
325 rtol, atol : float and array_like, optional
326 Relative and absolute tolerances. The solver keeps the local error
327 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
328 relative accuracy (number of correct digits), while `atol` controls
329 absolute accuracy (number of correct decimal places). To achieve the
330 desired `rtol`, set `atol` to be smaller than the smallest value that
331 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
332 allowable error. If `atol` is larger than ``rtol * abs(y)`` the
333 number of correct digits is not guaranteed. Conversely, to achieve the
334 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
335 than `atol`. If components of y have different scales, it might be
336 beneficial to set different `atol` values for different components by
337 passing array_like with shape (n,) for `atol`. Default values are
338 1e-3 for `rtol` and 1e-6 for `atol`.
339 vectorized : bool, optional
340 Whether `fun` is implemented in a vectorized fashion. Default is False.
342 Attributes
343 ----------
344 n : int
345 Number of equations.
346 status : string
347 Current status of the solver: 'running', 'finished' or 'failed'.
348 t_bound : float
349 Boundary time.
350 direction : float
351 Integration direction: +1 or -1.
352 t : float
353 Current time.
354 y : ndarray
355 Current state.
356 t_old : float
357 Previous time. None if no steps were made yet.
358 step_size : float
359 Size of the last successful step. None if no steps were made yet.
360 nfev : int
361 Number evaluations of the system's right-hand side.
362 njev : int
363 Number of evaluations of the Jacobian. Is always 0 for this solver as it does not use the Jacobian.
364 nlu : int
365 Number of LU decompositions. Is always 0 for this solver.
367 References
368 ----------
369 .. [1] J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
370 formulae", Journal of Computational and Applied Mathematics, Vol. 6,
371 No. 1, pp. 19-26, 1980.
372 .. [2] L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
373 of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
374 """
375 order = 5
376 error_estimator_order = 4
377 n_stages = 6
378 C = np.array([0, 1/5, 3/10, 4/5, 8/9, 1])
379 A = np.array([
380 [0, 0, 0, 0, 0],
381 [1/5, 0, 0, 0, 0],
382 [3/40, 9/40, 0, 0, 0],
383 [44/45, -56/15, 32/9, 0, 0],
384 [19372/6561, -25360/2187, 64448/6561, -212/729, 0],
385 [9017/3168, -355/33, 46732/5247, 49/176, -5103/18656]
386 ])
387 B = np.array([35/384, 0, 500/1113, 125/192, -2187/6784, 11/84])
388 E = np.array([-71/57600, 0, 71/16695, -71/1920, 17253/339200, -22/525,
389 1/40])
390 # Corresponds to the optimum value of c_6 from [2]_.
391 P = np.array([
392 [1, -8048581381/2820520608, 8663915743/2820520608,
393 -12715105075/11282082432],
394 [0, 0, 0, 0],
395 [0, 131558114200/32700410799, -68118460800/10900136933,
396 87487479700/32700410799],
397 [0, -1754552775/470086768, 14199869525/1410260304,
398 -10690763975/1880347072],
399 [0, 127303824393/49829197408, -318862633887/49829197408,
400 701980252875 / 199316789632],
401 [0, -282668133/205662961, 2019193451/616988883, -1453857185/822651844],
402 [0, 40617522/29380423, -110615467/29380423, 69997945/29380423]])
405class DOP853(RungeKutta):
406 """Explicit Runge-Kutta method of order 8.
408 This is a Python implementation of "DOP853" algorithm originally written
409 in Fortran [1]_, [2]_. Note that this is not a literate translation, but
410 the algorithmic core and coefficients are the same.
412 Can be applied in the complex domain.
414 Parameters
415 ----------
416 fun : callable
417 Right-hand side of the system. The calling signature is ``fun(t, y)``.
418 Here, ``t`` is a scalar, and there are two options for the ndarray ``y``:
419 It can either have shape (n,); then ``fun`` must return array_like with
420 shape (n,). Alternatively it can have shape (n, k); then ``fun``
421 must return an array_like with shape (n, k), i.e. each column
422 corresponds to a single column in ``y``. The choice between the two
423 options is determined by `vectorized` argument (see below).
424 t0 : float
425 Initial time.
426 y0 : array_like, shape (n,)
427 Initial state.
428 t_bound : float
429 Boundary time - the integration won't continue beyond it. It also
430 determines the direction of the integration.
431 first_step : float or None, optional
432 Initial step size. Default is ``None`` which means that the algorithm
433 should choose.
434 max_step : float, optional
435 Maximum allowed step size. Default is np.inf, i.e. the step size is not
436 bounded and determined solely by the solver.
437 rtol, atol : float and array_like, optional
438 Relative and absolute tolerances. The solver keeps the local error
439 estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
440 relative accuracy (number of correct digits), while `atol` controls
441 absolute accuracy (number of correct decimal places). To achieve the
442 desired `rtol`, set `atol` to be smaller than the smallest value that
443 can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
444 allowable error. If `atol` is larger than ``rtol * abs(y)`` the
445 number of correct digits is not guaranteed. Conversely, to achieve the
446 desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
447 than `atol`. If components of y have different scales, it might be
448 beneficial to set different `atol` values for different components by
449 passing array_like with shape (n,) for `atol`. Default values are
450 1e-3 for `rtol` and 1e-6 for `atol`.
451 vectorized : bool, optional
452 Whether `fun` is implemented in a vectorized fashion. Default is False.
454 Attributes
455 ----------
456 n : int
457 Number of equations.
458 status : string
459 Current status of the solver: 'running', 'finished' or 'failed'.
460 t_bound : float
461 Boundary time.
462 direction : float
463 Integration direction: +1 or -1.
464 t : float
465 Current time.
466 y : ndarray
467 Current state.
468 t_old : float
469 Previous time. None if no steps were made yet.
470 step_size : float
471 Size of the last successful step. None if no steps were made yet.
472 nfev : int
473 Number evaluations of the system's right-hand side.
474 njev : int
475 Number of evaluations of the Jacobian. Is always 0 for this solver
476 as it does not use the Jacobian.
477 nlu : int
478 Number of LU decompositions. Is always 0 for this solver.
480 References
481 ----------
482 .. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
483 Equations I: Nonstiff Problems", Sec. II.
484 .. [2] `Page with original Fortran code of DOP853
485 <http://www.unige.ch/~hairer/software.html>`_.
486 """
487 n_stages = dop853_coefficients.N_STAGES
488 order = 8
489 error_estimator_order = 7
490 A = dop853_coefficients.A[:n_stages, :n_stages]
491 B = dop853_coefficients.B
492 C = dop853_coefficients.C[:n_stages]
493 E3 = dop853_coefficients.E3
494 E5 = dop853_coefficients.E5
495 D = dop853_coefficients.D
497 A_EXTRA = dop853_coefficients.A[n_stages + 1:]
498 C_EXTRA = dop853_coefficients.C[n_stages + 1:]
500 def __init__(self, fun, t0, y0, t_bound, max_step=np.inf,
501 rtol=1e-3, atol=1e-6, vectorized=False,
502 first_step=None, **extraneous):
503 super().__init__(fun, t0, y0, t_bound, max_step, rtol, atol,
504 vectorized, first_step, **extraneous)
505 self.K_extended = np.empty((dop853_coefficients.N_STAGES_EXTENDED,
506 self.n), dtype=self.y.dtype)
507 self.K = self.K_extended[:self.n_stages + 1]
509 def _estimate_error(self, K, h): # Left for testing purposes.
510 err5 = np.dot(K.T, self.E5)
511 err3 = np.dot(K.T, self.E3)
512 denom = np.hypot(np.abs(err5), 0.1 * np.abs(err3))
513 correction_factor = np.ones_like(err5)
514 mask = denom > 0
515 correction_factor[mask] = np.abs(err5[mask]) / denom[mask]
516 return h * err5 * correction_factor
518 def _estimate_error_norm(self, K, h, scale):
519 err5 = np.dot(K.T, self.E5) / scale
520 err3 = np.dot(K.T, self.E3) / scale
521 err5_norm_2 = np.linalg.norm(err5)**2
522 err3_norm_2 = np.linalg.norm(err3)**2
523 if err5_norm_2 == 0 and err3_norm_2 == 0:
524 return 0.0
525 denom = err5_norm_2 + 0.01 * err3_norm_2
526 return np.abs(h) * err5_norm_2 / np.sqrt(denom * len(scale))
528 def _dense_output_impl(self):
529 K = self.K_extended
530 h = self.h_previous
531 for s, (a, c) in enumerate(zip(self.A_EXTRA, self.C_EXTRA),
532 start=self.n_stages + 1):
533 dy = np.dot(K[:s].T, a[:s]) * h
534 K[s] = self.fun(self.t_old + c * h, self.y_old + dy)
536 F = np.empty((dop853_coefficients.INTERPOLATOR_POWER, self.n),
537 dtype=self.y_old.dtype)
539 f_old = K[0]
540 delta_y = self.y - self.y_old
542 F[0] = delta_y
543 F[1] = h * f_old - delta_y
544 F[2] = 2 * delta_y - h * (self.f + f_old)
545 F[3:] = h * np.dot(self.D, K)
547 return Dop853DenseOutput(self.t_old, self.t, self.y_old, F)
550class RkDenseOutput(DenseOutput):
551 def __init__(self, t_old, t, y_old, Q):
552 super().__init__(t_old, t)
553 self.h = t - t_old
554 self.Q = Q
555 self.order = Q.shape[1] - 1
556 self.y_old = y_old
558 def _call_impl(self, t):
559 x = (t - self.t_old) / self.h
560 if t.ndim == 0:
561 p = np.tile(x, self.order + 1)
562 p = np.cumprod(p)
563 else:
564 p = np.tile(x, (self.order + 1, 1))
565 p = np.cumprod(p, axis=0)
566 y = self.h * np.dot(self.Q, p)
567 if y.ndim == 2:
568 y += self.y_old[:, None]
569 else:
570 y += self.y_old
572 return y
575class Dop853DenseOutput(DenseOutput):
576 def __init__(self, t_old, t, y_old, F):
577 super().__init__(t_old, t)
578 self.h = t - t_old
579 self.F = F
580 self.y_old = y_old
582 def _call_impl(self, t):
583 x = (t - self.t_old) / self.h
585 if t.ndim == 0:
586 y = np.zeros_like(self.y_old)
587 else:
588 x = x[:, None]
589 y = np.zeros((len(x), len(self.y_old)), dtype=self.y_old.dtype)
591 for i, f in enumerate(reversed(self.F)):
592 y += f
593 if i % 2 == 0:
594 y *= x
595 else:
596 y *= 1 - x
597 y += self.y_old
599 return y.T