Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_ivp/base.py: 18%
98 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
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 def fun_wrapped(t, y):
20 return np.asarray(fun(t, y), dtype=dtype)
22 return fun_wrapped, y0
25class OdeSolver:
26 """Base class for ODE solvers.
28 In order to implement a new solver you need to follow the guidelines:
30 1. A constructor must accept parameters presented in the base class
31 (listed below) along with any other parameters specific to a solver.
32 2. A constructor must accept arbitrary extraneous arguments
33 ``**extraneous``, but warn that these arguments are irrelevant
34 using `common.warn_extraneous` function. Do not pass these
35 arguments to the base class.
36 3. A solver must implement a private method `_step_impl(self)` which
37 propagates a solver one step further. It must return tuple
38 ``(success, message)``, where ``success`` is a boolean indicating
39 whether a step was successful, and ``message`` is a string
40 containing description of a failure if a step failed or None
41 otherwise.
42 4. A solver must implement a private method `_dense_output_impl(self)`,
43 which returns a `DenseOutput` object covering the last successful
44 step.
45 5. A solver must have attributes listed below in Attributes section.
46 Note that ``t_old`` and ``step_size`` are updated automatically.
47 6. Use `fun(self, t, y)` method for the system rhs evaluation, this
48 way the number of function evaluations (`nfev`) will be tracked
49 automatically.
50 7. For convenience, a base class provides `fun_single(self, t, y)` and
51 `fun_vectorized(self, t, y)` for evaluating the rhs in
52 non-vectorized and vectorized fashions respectively (regardless of
53 how `fun` from the constructor is implemented). These calls don't
54 increment `nfev`.
55 8. If a solver uses a Jacobian matrix and LU decompositions, it should
56 track the number of Jacobian evaluations (`njev`) and the number of
57 LU decompositions (`nlu`).
58 9. By convention, the function evaluations used to compute a finite
59 difference approximation of the Jacobian should not be counted in
60 `nfev`, thus use `fun_single(self, t, y)` or
61 `fun_vectorized(self, t, y)` when computing a finite difference
62 approximation of the Jacobian.
64 Parameters
65 ----------
66 fun : callable
67 Right-hand side of the system. The calling signature is ``fun(t, y)``.
68 Here ``t`` is a scalar and there are two options for ndarray ``y``.
69 It can either have shape (n,), then ``fun`` must return array_like with
70 shape (n,). Or, alternatively, it can have shape (n, n_points), then
71 ``fun`` must return array_like with shape (n, n_points) (each column
72 corresponds to a single column in ``y``). The choice between the two
73 options is determined by `vectorized` argument (see below).
74 t0 : float
75 Initial time.
76 y0 : array_like, shape (n,)
77 Initial state.
78 t_bound : float
79 Boundary time --- the integration won't continue beyond it. It also
80 determines the direction of the integration.
81 vectorized : bool
82 Whether `fun` is implemented in a vectorized fashion.
83 support_complex : bool, optional
84 Whether integration in a complex domain should be supported.
85 Generally determined by a derived solver class capabilities.
86 Default is False.
88 Attributes
89 ----------
90 n : int
91 Number of equations.
92 status : string
93 Current status of the solver: 'running', 'finished' or 'failed'.
94 t_bound : float
95 Boundary time.
96 direction : float
97 Integration direction: +1 or -1.
98 t : float
99 Current time.
100 y : ndarray
101 Current state.
102 t_old : float
103 Previous time. None if no steps were made yet.
104 step_size : float
105 Size of the last successful step. None if no steps were made yet.
106 nfev : int
107 Number of the system's rhs evaluations.
108 njev : int
109 Number of the Jacobian evaluations.
110 nlu : int
111 Number of LU decompositions.
112 """
113 TOO_SMALL_STEP = "Required step size is less than spacing between numbers."
115 def __init__(self, fun, t0, y0, t_bound, vectorized,
116 support_complex=False):
117 self.t_old = None
118 self.t = t0
119 self._fun, self.y = check_arguments(fun, y0, support_complex)
120 self.t_bound = t_bound
121 self.vectorized = vectorized
123 if vectorized:
124 def fun_single(t, y):
125 return self._fun(t, y[:, None]).ravel()
126 fun_vectorized = self._fun
127 else:
128 fun_single = self._fun
130 def fun_vectorized(t, y):
131 f = np.empty_like(y)
132 for i, yi in enumerate(y.T):
133 f[:, i] = self._fun(t, yi)
134 return f
136 def fun(t, y):
137 self.nfev += 1
138 return self.fun_single(t, y)
140 self.fun = fun
141 self.fun_single = fun_single
142 self.fun_vectorized = fun_vectorized
144 self.direction = np.sign(t_bound - t0) if t_bound != t0 else 1
145 self.n = self.y.size
146 self.status = 'running'
148 self.nfev = 0
149 self.njev = 0
150 self.nlu = 0
152 @property
153 def step_size(self):
154 if self.t_old is None:
155 return None
156 else:
157 return np.abs(self.t - self.t_old)
159 def step(self):
160 """Perform one integration step.
162 Returns
163 -------
164 message : string or None
165 Report from the solver. Typically a reason for a failure if
166 `self.status` is 'failed' after the step was taken or None
167 otherwise.
168 """
169 if self.status != 'running':
170 raise RuntimeError("Attempt to step on a failed or finished "
171 "solver.")
173 if self.n == 0 or self.t == self.t_bound:
174 # Handle corner cases of empty solver or no integration.
175 self.t_old = self.t
176 self.t = self.t_bound
177 message = None
178 self.status = 'finished'
179 else:
180 t = self.t
181 success, message = self._step_impl()
183 if not success:
184 self.status = 'failed'
185 else:
186 self.t_old = t
187 if self.direction * (self.t - self.t_bound) >= 0:
188 self.status = 'finished'
190 return message
192 def dense_output(self):
193 """Compute a local interpolant over the last successful step.
195 Returns
196 -------
197 sol : `DenseOutput`
198 Local interpolant over the last successful step.
199 """
200 if self.t_old is None:
201 raise RuntimeError("Dense output is available after a successful "
202 "step was made.")
204 if self.n == 0 or self.t == self.t_old:
205 # Handle corner cases of empty solver and no integration.
206 return ConstantDenseOutput(self.t_old, self.t, self.y)
207 else:
208 return self._dense_output_impl()
210 def _step_impl(self):
211 raise NotImplementedError
213 def _dense_output_impl(self):
214 raise NotImplementedError
217class DenseOutput:
218 """Base class for local interpolant over step made by an ODE solver.
220 It interpolates between `t_min` and `t_max` (see Attributes below).
221 Evaluation outside this interval is not forbidden, but the accuracy is not
222 guaranteed.
224 Attributes
225 ----------
226 t_min, t_max : float
227 Time range of the interpolation.
228 """
229 def __init__(self, t_old, t):
230 self.t_old = t_old
231 self.t = t
232 self.t_min = min(t, t_old)
233 self.t_max = max(t, t_old)
235 def __call__(self, t):
236 """Evaluate the interpolant.
238 Parameters
239 ----------
240 t : float or array_like with shape (n_points,)
241 Points to evaluate the solution at.
243 Returns
244 -------
245 y : ndarray, shape (n,) or (n, n_points)
246 Computed values. Shape depends on whether `t` was a scalar or a
247 1-D array.
248 """
249 t = np.asarray(t)
250 if t.ndim > 1:
251 raise ValueError("`t` must be a float or a 1-D array.")
252 return self._call_impl(t)
254 def _call_impl(self, t):
255 raise NotImplementedError
258class ConstantDenseOutput(DenseOutput):
259 """Constant value interpolator.
261 This class used for degenerate integration cases: equal integration limits
262 or a system with 0 equations.
263 """
264 def __init__(self, t_old, t, value):
265 super().__init__(t_old, t)
266 self.value = value
268 def _call_impl(self, t):
269 if t.ndim == 0:
270 return self.value
271 else:
272 ret = np.empty((self.value.shape[0], t.shape[0]))
273 ret[:] = self.value[:, None]
274 return ret