Coverage for /usr/lib/python3/dist-packages/scipy/integrate/_ivp/base.py: 18%
100 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
4def check_arguments(fun, y0, support_complex):
5 """Helper function for checking arguments common to all solvers."""
6 y0 = np.asarray(y0)
7 if np.issubdtype(y0.dtype, np.complexfloating):
8 if not support_complex:
9 raise ValueError("`y0` is complex, but the chosen solver does "
10 "not support integration in a complex domain.")
11 dtype = complex
12 else:
13 dtype = float
14 y0 = y0.astype(dtype, copy=False)
16 if y0.ndim != 1:
17 raise ValueError("`y0` must be 1-dimensional.")
19 if not np.isfinite(y0).all():
20 raise ValueError("All components of the initial state `y0` must be finite.")
22 def fun_wrapped(t, y):
23 return np.asarray(fun(t, y), dtype=dtype)
25 return fun_wrapped, y0
28class OdeSolver:
29 """Base class for ODE solvers.
31 In order to implement a new solver you need to follow the guidelines:
33 1. A constructor must accept parameters presented in the base class
34 (listed below) along with any other parameters specific to a solver.
35 2. A constructor must accept arbitrary extraneous arguments
36 ``**extraneous``, but warn that these arguments are irrelevant
37 using `common.warn_extraneous` function. Do not pass these
38 arguments to the base class.
39 3. A solver must implement a private method `_step_impl(self)` which
40 propagates a solver one step further. It must return tuple
41 ``(success, message)``, where ``success`` is a boolean indicating
42 whether a step was successful, and ``message`` is a string
43 containing description of a failure if a step failed or None
44 otherwise.
45 4. A solver must implement a private method `_dense_output_impl(self)`,
46 which returns a `DenseOutput` object covering the last successful
47 step.
48 5. A solver must have attributes listed below in Attributes section.
49 Note that ``t_old`` and ``step_size`` are updated automatically.
50 6. Use `fun(self, t, y)` method for the system rhs evaluation, this
51 way the number of function evaluations (`nfev`) will be tracked
52 automatically.
53 7. For convenience, a base class provides `fun_single(self, t, y)` and
54 `fun_vectorized(self, t, y)` for evaluating the rhs in
55 non-vectorized and vectorized fashions respectively (regardless of
56 how `fun` from the constructor is implemented). These calls don't
57 increment `nfev`.
58 8. If a solver uses a Jacobian matrix and LU decompositions, it should
59 track the number of Jacobian evaluations (`njev`) and the number of
60 LU decompositions (`nlu`).
61 9. By convention, the function evaluations used to compute a finite
62 difference approximation of the Jacobian should not be counted in
63 `nfev`, thus use `fun_single(self, t, y)` or
64 `fun_vectorized(self, t, y)` when computing a finite difference
65 approximation of the Jacobian.
67 Parameters
68 ----------
69 fun : callable
70 Right-hand side of the system: the time derivative of the state ``y``
71 at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a
72 scalar and ``y`` is an ndarray with ``len(y) = len(y0)``. ``fun`` must
73 return an array of the same shape as ``y``. See `vectorized` for more
74 information.
75 t0 : float
76 Initial time.
77 y0 : array_like, shape (n,)
78 Initial state.
79 t_bound : float
80 Boundary time --- the integration won't continue beyond it. It also
81 determines the direction of the integration.
82 vectorized : bool
83 Whether `fun` can be called in a vectorized fashion. Default is False.
85 If ``vectorized`` is False, `fun` will always be called with ``y`` of
86 shape ``(n,)``, where ``n = len(y0)``.
88 If ``vectorized`` is True, `fun` may be called with ``y`` of shape
89 ``(n, k)``, where ``k`` is an integer. In this case, `fun` must behave
90 such that ``fun(t, y)[:, i] == fun(t, y[:, i])`` (i.e. each column of
91 the returned array is the time derivative of the state corresponding
92 with a column of ``y``).
94 Setting ``vectorized=True`` allows for faster finite difference
95 approximation of the Jacobian by methods 'Radau' and 'BDF', but
96 will result in slower execution for other methods. It can also
97 result in slower overall execution for 'Radau' and 'BDF' in some
98 circumstances (e.g. small ``len(y0)``).
99 support_complex : bool, optional
100 Whether integration in a complex domain should be supported.
101 Generally determined by a derived solver class capabilities.
102 Default is False.
104 Attributes
105 ----------
106 n : int
107 Number of equations.
108 status : string
109 Current status of the solver: 'running', 'finished' or 'failed'.
110 t_bound : float
111 Boundary time.
112 direction : float
113 Integration direction: +1 or -1.
114 t : float
115 Current time.
116 y : ndarray
117 Current state.
118 t_old : float
119 Previous time. None if no steps were made yet.
120 step_size : float
121 Size of the last successful step. None if no steps were made yet.
122 nfev : int
123 Number of the system's rhs evaluations.
124 njev : int
125 Number of the Jacobian evaluations.
126 nlu : int
127 Number of LU decompositions.
128 """
129 TOO_SMALL_STEP = "Required step size is less than spacing between numbers."
131 def __init__(self, fun, t0, y0, t_bound, vectorized,
132 support_complex=False):
133 self.t_old = None
134 self.t = t0
135 self._fun, self.y = check_arguments(fun, y0, support_complex)
136 self.t_bound = t_bound
137 self.vectorized = vectorized
139 if vectorized:
140 def fun_single(t, y):
141 return self._fun(t, y[:, None]).ravel()
142 fun_vectorized = self._fun
143 else:
144 fun_single = self._fun
146 def fun_vectorized(t, y):
147 f = np.empty_like(y)
148 for i, yi in enumerate(y.T):
149 f[:, i] = self._fun(t, yi)
150 return f
152 def fun(t, y):
153 self.nfev += 1
154 return self.fun_single(t, y)
156 self.fun = fun
157 self.fun_single = fun_single
158 self.fun_vectorized = fun_vectorized
160 self.direction = np.sign(t_bound - t0) if t_bound != t0 else 1
161 self.n = self.y.size
162 self.status = 'running'
164 self.nfev = 0
165 self.njev = 0
166 self.nlu = 0
168 @property
169 def step_size(self):
170 if self.t_old is None:
171 return None
172 else:
173 return np.abs(self.t - self.t_old)
175 def step(self):
176 """Perform one integration step.
178 Returns
179 -------
180 message : string or None
181 Report from the solver. Typically a reason for a failure if
182 `self.status` is 'failed' after the step was taken or None
183 otherwise.
184 """
185 if self.status != 'running':
186 raise RuntimeError("Attempt to step on a failed or finished "
187 "solver.")
189 if self.n == 0 or self.t == self.t_bound:
190 # Handle corner cases of empty solver or no integration.
191 self.t_old = self.t
192 self.t = self.t_bound
193 message = None
194 self.status = 'finished'
195 else:
196 t = self.t
197 success, message = self._step_impl()
199 if not success:
200 self.status = 'failed'
201 else:
202 self.t_old = t
203 if self.direction * (self.t - self.t_bound) >= 0:
204 self.status = 'finished'
206 return message
208 def dense_output(self):
209 """Compute a local interpolant over the last successful step.
211 Returns
212 -------
213 sol : `DenseOutput`
214 Local interpolant over the last successful step.
215 """
216 if self.t_old is None:
217 raise RuntimeError("Dense output is available after a successful "
218 "step was made.")
220 if self.n == 0 or self.t == self.t_old:
221 # Handle corner cases of empty solver and no integration.
222 return ConstantDenseOutput(self.t_old, self.t, self.y)
223 else:
224 return self._dense_output_impl()
226 def _step_impl(self):
227 raise NotImplementedError
229 def _dense_output_impl(self):
230 raise NotImplementedError
233class DenseOutput:
234 """Base class for local interpolant over step made by an ODE solver.
236 It interpolates between `t_min` and `t_max` (see Attributes below).
237 Evaluation outside this interval is not forbidden, but the accuracy is not
238 guaranteed.
240 Attributes
241 ----------
242 t_min, t_max : float
243 Time range of the interpolation.
244 """
245 def __init__(self, t_old, t):
246 self.t_old = t_old
247 self.t = t
248 self.t_min = min(t, t_old)
249 self.t_max = max(t, t_old)
251 def __call__(self, t):
252 """Evaluate the interpolant.
254 Parameters
255 ----------
256 t : float or array_like with shape (n_points,)
257 Points to evaluate the solution at.
259 Returns
260 -------
261 y : ndarray, shape (n,) or (n, n_points)
262 Computed values. Shape depends on whether `t` was a scalar or a
263 1-D array.
264 """
265 t = np.asarray(t)
266 if t.ndim > 1:
267 raise ValueError("`t` must be a float or a 1-D array.")
268 return self._call_impl(t)
270 def _call_impl(self, t):
271 raise NotImplementedError
274class ConstantDenseOutput(DenseOutput):
275 """Constant value interpolator.
277 This class used for degenerate integration cases: equal integration limits
278 or a system with 0 equations.
279 """
280 def __init__(self, t_old, t, value):
281 super().__init__(t_old, t)
282 self.value = value
284 def _call_impl(self, t):
285 if t.ndim == 0:
286 return self.value
287 else:
288 ret = np.empty((self.value.shape[0], t.shape[0]))
289 ret[:] = self.value[:, None]
290 return ret