Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/integrate/_bvp.py: 7%
377 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Boundary value problem solver."""
2from warnings import warn
4import numpy as np
5from numpy.linalg import pinv
7from scipy.sparse import coo_matrix, csc_matrix
8from scipy.sparse.linalg import splu
9from scipy.optimize import OptimizeResult
12EPS = np.finfo(float).eps
15def estimate_fun_jac(fun, x, y, p, f0=None):
16 """Estimate derivatives of an ODE system rhs with forward differences.
18 Returns
19 -------
20 df_dy : ndarray, shape (n, n, m)
21 Derivatives with respect to y. An element (i, j, q) corresponds to
22 d f_i(x_q, y_q) / d (y_q)_j.
23 df_dp : ndarray with shape (n, k, m) or None
24 Derivatives with respect to p. An element (i, j, q) corresponds to
25 d f_i(x_q, y_q, p) / d p_j. If `p` is empty, None is returned.
26 """
27 n, m = y.shape
28 if f0 is None:
29 f0 = fun(x, y, p)
31 dtype = y.dtype
33 df_dy = np.empty((n, n, m), dtype=dtype)
34 h = EPS**0.5 * (1 + np.abs(y))
35 for i in range(n):
36 y_new = y.copy()
37 y_new[i] += h[i]
38 hi = y_new[i] - y[i]
39 f_new = fun(x, y_new, p)
40 df_dy[:, i, :] = (f_new - f0) / hi
42 k = p.shape[0]
43 if k == 0:
44 df_dp = None
45 else:
46 df_dp = np.empty((n, k, m), dtype=dtype)
47 h = EPS**0.5 * (1 + np.abs(p))
48 for i in range(k):
49 p_new = p.copy()
50 p_new[i] += h[i]
51 hi = p_new[i] - p[i]
52 f_new = fun(x, y, p_new)
53 df_dp[:, i, :] = (f_new - f0) / hi
55 return df_dy, df_dp
58def estimate_bc_jac(bc, ya, yb, p, bc0=None):
59 """Estimate derivatives of boundary conditions with forward differences.
61 Returns
62 -------
63 dbc_dya : ndarray, shape (n + k, n)
64 Derivatives with respect to ya. An element (i, j) corresponds to
65 d bc_i / d ya_j.
66 dbc_dyb : ndarray, shape (n + k, n)
67 Derivatives with respect to yb. An element (i, j) corresponds to
68 d bc_i / d ya_j.
69 dbc_dp : ndarray with shape (n + k, k) or None
70 Derivatives with respect to p. An element (i, j) corresponds to
71 d bc_i / d p_j. If `p` is empty, None is returned.
72 """
73 n = ya.shape[0]
74 k = p.shape[0]
76 if bc0 is None:
77 bc0 = bc(ya, yb, p)
79 dtype = ya.dtype
81 dbc_dya = np.empty((n, n + k), dtype=dtype)
82 h = EPS**0.5 * (1 + np.abs(ya))
83 for i in range(n):
84 ya_new = ya.copy()
85 ya_new[i] += h[i]
86 hi = ya_new[i] - ya[i]
87 bc_new = bc(ya_new, yb, p)
88 dbc_dya[i] = (bc_new - bc0) / hi
89 dbc_dya = dbc_dya.T
91 h = EPS**0.5 * (1 + np.abs(yb))
92 dbc_dyb = np.empty((n, n + k), dtype=dtype)
93 for i in range(n):
94 yb_new = yb.copy()
95 yb_new[i] += h[i]
96 hi = yb_new[i] - yb[i]
97 bc_new = bc(ya, yb_new, p)
98 dbc_dyb[i] = (bc_new - bc0) / hi
99 dbc_dyb = dbc_dyb.T
101 if k == 0:
102 dbc_dp = None
103 else:
104 h = EPS**0.5 * (1 + np.abs(p))
105 dbc_dp = np.empty((k, n + k), dtype=dtype)
106 for i in range(k):
107 p_new = p.copy()
108 p_new[i] += h[i]
109 hi = p_new[i] - p[i]
110 bc_new = bc(ya, yb, p_new)
111 dbc_dp[i] = (bc_new - bc0) / hi
112 dbc_dp = dbc_dp.T
114 return dbc_dya, dbc_dyb, dbc_dp
117def compute_jac_indices(n, m, k):
118 """Compute indices for the collocation system Jacobian construction.
120 See `construct_global_jac` for the explanation.
121 """
122 i_col = np.repeat(np.arange((m - 1) * n), n)
123 j_col = (np.tile(np.arange(n), n * (m - 1)) +
124 np.repeat(np.arange(m - 1) * n, n**2))
126 i_bc = np.repeat(np.arange((m - 1) * n, m * n + k), n)
127 j_bc = np.tile(np.arange(n), n + k)
129 i_p_col = np.repeat(np.arange((m - 1) * n), k)
130 j_p_col = np.tile(np.arange(m * n, m * n + k), (m - 1) * n)
132 i_p_bc = np.repeat(np.arange((m - 1) * n, m * n + k), k)
133 j_p_bc = np.tile(np.arange(m * n, m * n + k), n + k)
135 i = np.hstack((i_col, i_col, i_bc, i_bc, i_p_col, i_p_bc))
136 j = np.hstack((j_col, j_col + n,
137 j_bc, j_bc + (m - 1) * n,
138 j_p_col, j_p_bc))
140 return i, j
143def stacked_matmul(a, b):
144 """Stacked matrix multiply: out[i,:,:] = np.dot(a[i,:,:], b[i,:,:]).
146 Empirical optimization. Use outer Python loop and BLAS for large
147 matrices, otherwise use a single einsum call.
148 """
149 if a.shape[1] > 50:
150 out = np.empty((a.shape[0], a.shape[1], b.shape[2]))
151 for i in range(a.shape[0]):
152 out[i] = np.dot(a[i], b[i])
153 return out
154 else:
155 return np.einsum('...ij,...jk->...ik', a, b)
158def construct_global_jac(n, m, k, i_jac, j_jac, h, df_dy, df_dy_middle, df_dp,
159 df_dp_middle, dbc_dya, dbc_dyb, dbc_dp):
160 """Construct the Jacobian of the collocation system.
162 There are n * m + k functions: m - 1 collocations residuals, each
163 containing n components, followed by n + k boundary condition residuals.
165 There are n * m + k variables: m vectors of y, each containing n
166 components, followed by k values of vector p.
168 For example, let m = 4, n = 2 and k = 1, then the Jacobian will have
169 the following sparsity structure:
171 1 1 2 2 0 0 0 0 5
172 1 1 2 2 0 0 0 0 5
173 0 0 1 1 2 2 0 0 5
174 0 0 1 1 2 2 0 0 5
175 0 0 0 0 1 1 2 2 5
176 0 0 0 0 1 1 2 2 5
178 3 3 0 0 0 0 4 4 6
179 3 3 0 0 0 0 4 4 6
180 3 3 0 0 0 0 4 4 6
182 Zeros denote identically zero values, other values denote different kinds
183 of blocks in the matrix (see below). The blank row indicates the separation
184 of collocation residuals from boundary conditions. And the blank column
185 indicates the separation of y values from p values.
187 Refer to [1]_ (p. 306) for the formula of n x n blocks for derivatives
188 of collocation residuals with respect to y.
190 Parameters
191 ----------
192 n : int
193 Number of equations in the ODE system.
194 m : int
195 Number of nodes in the mesh.
196 k : int
197 Number of the unknown parameters.
198 i_jac, j_jac : ndarray
199 Row and column indices returned by `compute_jac_indices`. They
200 represent different blocks in the Jacobian matrix in the following
201 order (see the scheme above):
203 * 1: m - 1 diagonal n x n blocks for the collocation residuals.
204 * 2: m - 1 off-diagonal n x n blocks for the collocation residuals.
205 * 3 : (n + k) x n block for the dependency of the boundary
206 conditions on ya.
207 * 4: (n + k) x n block for the dependency of the boundary
208 conditions on yb.
209 * 5: (m - 1) * n x k block for the dependency of the collocation
210 residuals on p.
211 * 6: (n + k) x k block for the dependency of the boundary
212 conditions on p.
214 df_dy : ndarray, shape (n, n, m)
215 Jacobian of f with respect to y computed at the mesh nodes.
216 df_dy_middle : ndarray, shape (n, n, m - 1)
217 Jacobian of f with respect to y computed at the middle between the
218 mesh nodes.
219 df_dp : ndarray with shape (n, k, m) or None
220 Jacobian of f with respect to p computed at the mesh nodes.
221 df_dp_middle : ndarray with shape (n, k, m - 1) or None
222 Jacobian of f with respect to p computed at the middle between the
223 mesh nodes.
224 dbc_dya, dbc_dyb : ndarray, shape (n, n)
225 Jacobian of bc with respect to ya and yb.
226 dbc_dp : ndarray with shape (n, k) or None
227 Jacobian of bc with respect to p.
229 Returns
230 -------
231 J : csc_matrix, shape (n * m + k, n * m + k)
232 Jacobian of the collocation system in a sparse form.
234 References
235 ----------
236 .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
237 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
238 Number 3, pp. 299-316, 2001.
239 """
240 df_dy = np.transpose(df_dy, (2, 0, 1))
241 df_dy_middle = np.transpose(df_dy_middle, (2, 0, 1))
243 h = h[:, np.newaxis, np.newaxis]
245 dtype = df_dy.dtype
247 # Computing diagonal n x n blocks.
248 dPhi_dy_0 = np.empty((m - 1, n, n), dtype=dtype)
249 dPhi_dy_0[:] = -np.identity(n)
250 dPhi_dy_0 -= h / 6 * (df_dy[:-1] + 2 * df_dy_middle)
251 T = stacked_matmul(df_dy_middle, df_dy[:-1])
252 dPhi_dy_0 -= h**2 / 12 * T
254 # Computing off-diagonal n x n blocks.
255 dPhi_dy_1 = np.empty((m - 1, n, n), dtype=dtype)
256 dPhi_dy_1[:] = np.identity(n)
257 dPhi_dy_1 -= h / 6 * (df_dy[1:] + 2 * df_dy_middle)
258 T = stacked_matmul(df_dy_middle, df_dy[1:])
259 dPhi_dy_1 += h**2 / 12 * T
261 values = np.hstack((dPhi_dy_0.ravel(), dPhi_dy_1.ravel(), dbc_dya.ravel(),
262 dbc_dyb.ravel()))
264 if k > 0:
265 df_dp = np.transpose(df_dp, (2, 0, 1))
266 df_dp_middle = np.transpose(df_dp_middle, (2, 0, 1))
267 T = stacked_matmul(df_dy_middle, df_dp[:-1] - df_dp[1:])
268 df_dp_middle += 0.125 * h * T
269 dPhi_dp = -h/6 * (df_dp[:-1] + df_dp[1:] + 4 * df_dp_middle)
270 values = np.hstack((values, dPhi_dp.ravel(), dbc_dp.ravel()))
272 J = coo_matrix((values, (i_jac, j_jac)))
273 return csc_matrix(J)
276def collocation_fun(fun, y, p, x, h):
277 """Evaluate collocation residuals.
279 This function lies in the core of the method. The solution is sought
280 as a cubic C1 continuous spline with derivatives matching the ODE rhs
281 at given nodes `x`. Collocation conditions are formed from the equality
282 of the spline derivatives and rhs of the ODE system in the middle points
283 between nodes.
285 Such method is classified to Lobbato IIIA family in ODE literature.
286 Refer to [1]_ for the formula and some discussion.
288 Returns
289 -------
290 col_res : ndarray, shape (n, m - 1)
291 Collocation residuals at the middle points of the mesh intervals.
292 y_middle : ndarray, shape (n, m - 1)
293 Values of the cubic spline evaluated at the middle points of the mesh
294 intervals.
295 f : ndarray, shape (n, m)
296 RHS of the ODE system evaluated at the mesh nodes.
297 f_middle : ndarray, shape (n, m - 1)
298 RHS of the ODE system evaluated at the middle points of the mesh
299 intervals (and using `y_middle`).
301 References
302 ----------
303 .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
304 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
305 Number 3, pp. 299-316, 2001.
306 """
307 f = fun(x, y, p)
308 y_middle = (0.5 * (y[:, 1:] + y[:, :-1]) -
309 0.125 * h * (f[:, 1:] - f[:, :-1]))
310 f_middle = fun(x[:-1] + 0.5 * h, y_middle, p)
311 col_res = y[:, 1:] - y[:, :-1] - h / 6 * (f[:, :-1] + f[:, 1:] +
312 4 * f_middle)
314 return col_res, y_middle, f, f_middle
317def prepare_sys(n, m, k, fun, bc, fun_jac, bc_jac, x, h):
318 """Create the function and the Jacobian for the collocation system."""
319 x_middle = x[:-1] + 0.5 * h
320 i_jac, j_jac = compute_jac_indices(n, m, k)
322 def col_fun(y, p):
323 return collocation_fun(fun, y, p, x, h)
325 def sys_jac(y, p, y_middle, f, f_middle, bc0):
326 if fun_jac is None:
327 df_dy, df_dp = estimate_fun_jac(fun, x, y, p, f)
328 df_dy_middle, df_dp_middle = estimate_fun_jac(
329 fun, x_middle, y_middle, p, f_middle)
330 else:
331 df_dy, df_dp = fun_jac(x, y, p)
332 df_dy_middle, df_dp_middle = fun_jac(x_middle, y_middle, p)
334 if bc_jac is None:
335 dbc_dya, dbc_dyb, dbc_dp = estimate_bc_jac(bc, y[:, 0], y[:, -1],
336 p, bc0)
337 else:
338 dbc_dya, dbc_dyb, dbc_dp = bc_jac(y[:, 0], y[:, -1], p)
340 return construct_global_jac(n, m, k, i_jac, j_jac, h, df_dy,
341 df_dy_middle, df_dp, df_dp_middle, dbc_dya,
342 dbc_dyb, dbc_dp)
344 return col_fun, sys_jac
347def solve_newton(n, m, h, col_fun, bc, jac, y, p, B, bvp_tol, bc_tol):
348 """Solve the nonlinear collocation system by a Newton method.
350 This is a simple Newton method with a backtracking line search. As
351 advised in [1]_, an affine-invariant criterion function F = ||J^-1 r||^2
352 is used, where J is the Jacobian matrix at the current iteration and r is
353 the vector or collocation residuals (values of the system lhs).
355 The method alters between full Newton iterations and the fixed-Jacobian
356 iterations based
358 There are other tricks proposed in [1]_, but they are not used as they
359 don't seem to improve anything significantly, and even break the
360 convergence on some test problems I tried.
362 All important parameters of the algorithm are defined inside the function.
364 Parameters
365 ----------
366 n : int
367 Number of equations in the ODE system.
368 m : int
369 Number of nodes in the mesh.
370 h : ndarray, shape (m-1,)
371 Mesh intervals.
372 col_fun : callable
373 Function computing collocation residuals.
374 bc : callable
375 Function computing boundary condition residuals.
376 jac : callable
377 Function computing the Jacobian of the whole system (including
378 collocation and boundary condition residuals). It is supposed to
379 return csc_matrix.
380 y : ndarray, shape (n, m)
381 Initial guess for the function values at the mesh nodes.
382 p : ndarray, shape (k,)
383 Initial guess for the unknown parameters.
384 B : ndarray with shape (n, n) or None
385 Matrix to force the S y(a) = 0 condition for a problems with the
386 singular term. If None, the singular term is assumed to be absent.
387 bvp_tol : float
388 Tolerance to which we want to solve a BVP.
389 bc_tol : float
390 Tolerance to which we want to satisfy the boundary conditions.
392 Returns
393 -------
394 y : ndarray, shape (n, m)
395 Final iterate for the function values at the mesh nodes.
396 p : ndarray, shape (k,)
397 Final iterate for the unknown parameters.
398 singular : bool
399 True, if the LU decomposition failed because Jacobian turned out
400 to be singular.
402 References
403 ----------
404 .. [1] U. Ascher, R. Mattheij and R. Russell "Numerical Solution of
405 Boundary Value Problems for Ordinary Differential Equations"
406 """
407 # We know that the solution residuals at the middle points of the mesh
408 # are connected with collocation residuals r_middle = 1.5 * col_res / h.
409 # As our BVP solver tries to decrease relative residuals below a certain
410 # tolerance, it seems reasonable to terminated Newton iterations by
411 # comparison of r_middle / (1 + np.abs(f_middle)) with a certain threshold,
412 # which we choose to be 1.5 orders lower than the BVP tolerance. We rewrite
413 # the condition as col_res < tol_r * (1 + np.abs(f_middle)), then tol_r
414 # should be computed as follows:
415 tol_r = 2/3 * h * 5e-2 * bvp_tol
417 # Maximum allowed number of Jacobian evaluation and factorization, in
418 # other words, the maximum number of full Newton iterations. A small value
419 # is recommended in the literature.
420 max_njev = 4
422 # Maximum number of iterations, considering that some of them can be
423 # performed with the fixed Jacobian. In theory, such iterations are cheap,
424 # but it's not that simple in Python.
425 max_iter = 8
427 # Minimum relative improvement of the criterion function to accept the
428 # step (Armijo constant).
429 sigma = 0.2
431 # Step size decrease factor for backtracking.
432 tau = 0.5
434 # Maximum number of backtracking steps, the minimum step is then
435 # tau ** n_trial.
436 n_trial = 4
438 col_res, y_middle, f, f_middle = col_fun(y, p)
439 bc_res = bc(y[:, 0], y[:, -1], p)
440 res = np.hstack((col_res.ravel(order='F'), bc_res))
442 njev = 0
443 singular = False
444 recompute_jac = True
445 for iteration in range(max_iter):
446 if recompute_jac:
447 J = jac(y, p, y_middle, f, f_middle, bc_res)
448 njev += 1
449 try:
450 LU = splu(J)
451 except RuntimeError:
452 singular = True
453 break
455 step = LU.solve(res)
456 cost = np.dot(step, step)
458 y_step = step[:m * n].reshape((n, m), order='F')
459 p_step = step[m * n:]
461 alpha = 1
462 for trial in range(n_trial + 1):
463 y_new = y - alpha * y_step
464 if B is not None:
465 y_new[:, 0] = np.dot(B, y_new[:, 0])
466 p_new = p - alpha * p_step
468 col_res, y_middle, f, f_middle = col_fun(y_new, p_new)
469 bc_res = bc(y_new[:, 0], y_new[:, -1], p_new)
470 res = np.hstack((col_res.ravel(order='F'), bc_res))
472 step_new = LU.solve(res)
473 cost_new = np.dot(step_new, step_new)
474 if cost_new < (1 - 2 * alpha * sigma) * cost:
475 break
477 if trial < n_trial:
478 alpha *= tau
480 y = y_new
481 p = p_new
483 if njev == max_njev:
484 break
486 if (np.all(np.abs(col_res) < tol_r * (1 + np.abs(f_middle))) and
487 np.all(np.abs(bc_res) < bc_tol)):
488 break
490 # If the full step was taken, then we are going to continue with
491 # the same Jacobian. This is the approach of BVP_SOLVER.
492 if alpha == 1:
493 step = step_new
494 cost = cost_new
495 recompute_jac = False
496 else:
497 recompute_jac = True
499 return y, p, singular
502def print_iteration_header():
503 print("{:^15}{:^15}{:^15}{:^15}{:^15}".format(
504 "Iteration", "Max residual", "Max BC residual", "Total nodes",
505 "Nodes added"))
508def print_iteration_progress(iteration, residual, bc_residual, total_nodes,
509 nodes_added):
510 print("{:^15}{:^15.2e}{:^15.2e}{:^15}{:^15}".format(
511 iteration, residual, bc_residual, total_nodes, nodes_added))
514class BVPResult(OptimizeResult):
515 pass
518TERMINATION_MESSAGES = {
519 0: "The algorithm converged to the desired accuracy.",
520 1: "The maximum number of mesh nodes is exceeded.",
521 2: "A singular Jacobian encountered when solving the collocation system.",
522 3: "The solver was unable to satisfy boundary conditions tolerance on iteration 10."
523}
526def estimate_rms_residuals(fun, sol, x, h, p, r_middle, f_middle):
527 """Estimate rms values of collocation residuals using Lobatto quadrature.
529 The residuals are defined as the difference between the derivatives of
530 our solution and rhs of the ODE system. We use relative residuals, i.e.,
531 normalized by 1 + np.abs(f). RMS values are computed as sqrt from the
532 normalized integrals of the squared relative residuals over each interval.
533 Integrals are estimated using 5-point Lobatto quadrature [1]_, we use the
534 fact that residuals at the mesh nodes are identically zero.
536 In [2] they don't normalize integrals by interval lengths, which gives
537 a higher rate of convergence of the residuals by the factor of h**0.5.
538 I chose to do such normalization for an ease of interpretation of return
539 values as RMS estimates.
541 Returns
542 -------
543 rms_res : ndarray, shape (m - 1,)
544 Estimated rms values of the relative residuals over each interval.
546 References
547 ----------
548 .. [1] http://mathworld.wolfram.com/LobattoQuadrature.html
549 .. [2] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
550 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
551 Number 3, pp. 299-316, 2001.
552 """
553 x_middle = x[:-1] + 0.5 * h
554 s = 0.5 * h * (3/7)**0.5
555 x1 = x_middle + s
556 x2 = x_middle - s
557 y1 = sol(x1)
558 y2 = sol(x2)
559 y1_prime = sol(x1, 1)
560 y2_prime = sol(x2, 1)
561 f1 = fun(x1, y1, p)
562 f2 = fun(x2, y2, p)
563 r1 = y1_prime - f1
564 r2 = y2_prime - f2
566 r_middle /= 1 + np.abs(f_middle)
567 r1 /= 1 + np.abs(f1)
568 r2 /= 1 + np.abs(f2)
570 r1 = np.sum(np.real(r1 * np.conj(r1)), axis=0)
571 r2 = np.sum(np.real(r2 * np.conj(r2)), axis=0)
572 r_middle = np.sum(np.real(r_middle * np.conj(r_middle)), axis=0)
574 return (0.5 * (32 / 45 * r_middle + 49 / 90 * (r1 + r2))) ** 0.5
577def create_spline(y, yp, x, h):
578 """Create a cubic spline given values and derivatives.
580 Formulas for the coefficients are taken from interpolate.CubicSpline.
582 Returns
583 -------
584 sol : PPoly
585 Constructed spline as a PPoly instance.
586 """
587 from scipy.interpolate import PPoly
589 n, m = y.shape
590 c = np.empty((4, n, m - 1), dtype=y.dtype)
591 slope = (y[:, 1:] - y[:, :-1]) / h
592 t = (yp[:, :-1] + yp[:, 1:] - 2 * slope) / h
593 c[0] = t / h
594 c[1] = (slope - yp[:, :-1]) / h - t
595 c[2] = yp[:, :-1]
596 c[3] = y[:, :-1]
597 c = np.moveaxis(c, 1, 0)
599 return PPoly(c, x, extrapolate=True, axis=1)
602def modify_mesh(x, insert_1, insert_2):
603 """Insert nodes into a mesh.
605 Nodes removal logic is not established, its impact on the solver is
606 presumably negligible. So, only insertion is done in this function.
608 Parameters
609 ----------
610 x : ndarray, shape (m,)
611 Mesh nodes.
612 insert_1 : ndarray
613 Intervals to each insert 1 new node in the middle.
614 insert_2 : ndarray
615 Intervals to each insert 2 new nodes, such that divide an interval
616 into 3 equal parts.
618 Returns
619 -------
620 x_new : ndarray
621 New mesh nodes.
623 Notes
624 -----
625 `insert_1` and `insert_2` should not have common values.
626 """
627 # Because np.insert implementation apparently varies with a version of
628 # NumPy, we use a simple and reliable approach with sorting.
629 return np.sort(np.hstack((
630 x,
631 0.5 * (x[insert_1] + x[insert_1 + 1]),
632 (2 * x[insert_2] + x[insert_2 + 1]) / 3,
633 (x[insert_2] + 2 * x[insert_2 + 1]) / 3
634 )))
637def wrap_functions(fun, bc, fun_jac, bc_jac, k, a, S, D, dtype):
638 """Wrap functions for unified usage in the solver."""
639 if fun_jac is None:
640 fun_jac_wrapped = None
642 if bc_jac is None:
643 bc_jac_wrapped = None
645 if k == 0:
646 def fun_p(x, y, _):
647 return np.asarray(fun(x, y), dtype)
649 def bc_wrapped(ya, yb, _):
650 return np.asarray(bc(ya, yb), dtype)
652 if fun_jac is not None:
653 def fun_jac_p(x, y, _):
654 return np.asarray(fun_jac(x, y), dtype), None
656 if bc_jac is not None:
657 def bc_jac_wrapped(ya, yb, _):
658 dbc_dya, dbc_dyb = bc_jac(ya, yb)
659 return (np.asarray(dbc_dya, dtype),
660 np.asarray(dbc_dyb, dtype), None)
661 else:
662 def fun_p(x, y, p):
663 return np.asarray(fun(x, y, p), dtype)
665 def bc_wrapped(x, y, p):
666 return np.asarray(bc(x, y, p), dtype)
668 if fun_jac is not None:
669 def fun_jac_p(x, y, p):
670 df_dy, df_dp = fun_jac(x, y, p)
671 return np.asarray(df_dy, dtype), np.asarray(df_dp, dtype)
673 if bc_jac is not None:
674 def bc_jac_wrapped(ya, yb, p):
675 dbc_dya, dbc_dyb, dbc_dp = bc_jac(ya, yb, p)
676 return (np.asarray(dbc_dya, dtype), np.asarray(dbc_dyb, dtype),
677 np.asarray(dbc_dp, dtype))
679 if S is None:
680 fun_wrapped = fun_p
681 else:
682 def fun_wrapped(x, y, p):
683 f = fun_p(x, y, p)
684 if x[0] == a:
685 f[:, 0] = np.dot(D, f[:, 0])
686 f[:, 1:] += np.dot(S, y[:, 1:]) / (x[1:] - a)
687 else:
688 f += np.dot(S, y) / (x - a)
689 return f
691 if fun_jac is not None:
692 if S is None:
693 fun_jac_wrapped = fun_jac_p
694 else:
695 Sr = S[:, :, np.newaxis]
697 def fun_jac_wrapped(x, y, p):
698 df_dy, df_dp = fun_jac_p(x, y, p)
699 if x[0] == a:
700 df_dy[:, :, 0] = np.dot(D, df_dy[:, :, 0])
701 df_dy[:, :, 1:] += Sr / (x[1:] - a)
702 else:
703 df_dy += Sr / (x - a)
705 return df_dy, df_dp
707 return fun_wrapped, bc_wrapped, fun_jac_wrapped, bc_jac_wrapped
710def solve_bvp(fun, bc, x, y, p=None, S=None, fun_jac=None, bc_jac=None,
711 tol=1e-3, max_nodes=1000, verbose=0, bc_tol=None):
712 """Solve a boundary value problem for a system of ODEs.
714 This function numerically solves a first order system of ODEs subject to
715 two-point boundary conditions::
717 dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b
718 bc(y(a), y(b), p) = 0
720 Here x is a 1-D independent variable, y(x) is an N-D
721 vector-valued function and p is a k-D vector of unknown
722 parameters which is to be found along with y(x). For the problem to be
723 determined, there must be n + k boundary conditions, i.e., bc must be an
724 (n + k)-D function.
726 The last singular term on the right-hand side of the system is optional.
727 It is defined by an n-by-n matrix S, such that the solution must satisfy
728 S y(a) = 0. This condition will be forced during iterations, so it must not
729 contradict boundary conditions. See [2]_ for the explanation how this term
730 is handled when solving BVPs numerically.
732 Problems in a complex domain can be solved as well. In this case, y and p
733 are considered to be complex, and f and bc are assumed to be complex-valued
734 functions, but x stays real. Note that f and bc must be complex
735 differentiable (satisfy Cauchy-Riemann equations [4]_), otherwise you
736 should rewrite your problem for real and imaginary parts separately. To
737 solve a problem in a complex domain, pass an initial guess for y with a
738 complex data type (see below).
740 Parameters
741 ----------
742 fun : callable
743 Right-hand side of the system. The calling signature is ``fun(x, y)``,
744 or ``fun(x, y, p)`` if parameters are present. All arguments are
745 ndarray: ``x`` with shape (m,), ``y`` with shape (n, m), meaning that
746 ``y[:, i]`` corresponds to ``x[i]``, and ``p`` with shape (k,). The
747 return value must be an array with shape (n, m) and with the same
748 layout as ``y``.
749 bc : callable
750 Function evaluating residuals of the boundary conditions. The calling
751 signature is ``bc(ya, yb)``, or ``bc(ya, yb, p)`` if parameters are
752 present. All arguments are ndarray: ``ya`` and ``yb`` with shape (n,),
753 and ``p`` with shape (k,). The return value must be an array with
754 shape (n + k,).
755 x : array_like, shape (m,)
756 Initial mesh. Must be a strictly increasing sequence of real numbers
757 with ``x[0]=a`` and ``x[-1]=b``.
758 y : array_like, shape (n, m)
759 Initial guess for the function values at the mesh nodes, ith column
760 corresponds to ``x[i]``. For problems in a complex domain pass `y`
761 with a complex data type (even if the initial guess is purely real).
762 p : array_like with shape (k,) or None, optional
763 Initial guess for the unknown parameters. If None (default), it is
764 assumed that the problem doesn't depend on any parameters.
765 S : array_like with shape (n, n) or None
766 Matrix defining the singular term. If None (default), the problem is
767 solved without the singular term.
768 fun_jac : callable or None, optional
769 Function computing derivatives of f with respect to y and p. The
770 calling signature is ``fun_jac(x, y)``, or ``fun_jac(x, y, p)`` if
771 parameters are present. The return must contain 1 or 2 elements in the
772 following order:
774 * df_dy : array_like with shape (n, n, m), where an element
775 (i, j, q) equals to d f_i(x_q, y_q, p) / d (y_q)_j.
776 * df_dp : array_like with shape (n, k, m), where an element
777 (i, j, q) equals to d f_i(x_q, y_q, p) / d p_j.
779 Here q numbers nodes at which x and y are defined, whereas i and j
780 number vector components. If the problem is solved without unknown
781 parameters, df_dp should not be returned.
783 If `fun_jac` is None (default), the derivatives will be estimated
784 by the forward finite differences.
785 bc_jac : callable or None, optional
786 Function computing derivatives of bc with respect to ya, yb, and p.
787 The calling signature is ``bc_jac(ya, yb)``, or ``bc_jac(ya, yb, p)``
788 if parameters are present. The return must contain 2 or 3 elements in
789 the following order:
791 * dbc_dya : array_like with shape (n, n), where an element (i, j)
792 equals to d bc_i(ya, yb, p) / d ya_j.
793 * dbc_dyb : array_like with shape (n, n), where an element (i, j)
794 equals to d bc_i(ya, yb, p) / d yb_j.
795 * dbc_dp : array_like with shape (n, k), where an element (i, j)
796 equals to d bc_i(ya, yb, p) / d p_j.
798 If the problem is solved without unknown parameters, dbc_dp should not
799 be returned.
801 If `bc_jac` is None (default), the derivatives will be estimated by
802 the forward finite differences.
803 tol : float, optional
804 Desired tolerance of the solution. If we define ``r = y' - f(x, y)``,
805 where y is the found solution, then the solver tries to achieve on each
806 mesh interval ``norm(r / (1 + abs(f)) < tol``, where ``norm`` is
807 estimated in a root mean squared sense (using a numerical quadrature
808 formula). Default is 1e-3.
809 max_nodes : int, optional
810 Maximum allowed number of the mesh nodes. If exceeded, the algorithm
811 terminates. Default is 1000.
812 verbose : {0, 1, 2}, optional
813 Level of algorithm's verbosity:
815 * 0 (default) : work silently.
816 * 1 : display a termination report.
817 * 2 : display progress during iterations.
818 bc_tol : float, optional
819 Desired absolute tolerance for the boundary condition residuals: `bc`
820 value should satisfy ``abs(bc) < bc_tol`` component-wise.
821 Equals to `tol` by default. Up to 10 iterations are allowed to achieve this
822 tolerance.
824 Returns
825 -------
826 Bunch object with the following fields defined:
827 sol : PPoly
828 Found solution for y as `scipy.interpolate.PPoly` instance, a C1
829 continuous cubic spline.
830 p : ndarray or None, shape (k,)
831 Found parameters. None, if the parameters were not present in the
832 problem.
833 x : ndarray, shape (m,)
834 Nodes of the final mesh.
835 y : ndarray, shape (n, m)
836 Solution values at the mesh nodes.
837 yp : ndarray, shape (n, m)
838 Solution derivatives at the mesh nodes.
839 rms_residuals : ndarray, shape (m - 1,)
840 RMS values of the relative residuals over each mesh interval (see the
841 description of `tol` parameter).
842 niter : int
843 Number of completed iterations.
844 status : int
845 Reason for algorithm termination:
847 * 0: The algorithm converged to the desired accuracy.
848 * 1: The maximum number of mesh nodes is exceeded.
849 * 2: A singular Jacobian encountered when solving the collocation
850 system.
852 message : string
853 Verbal description of the termination reason.
854 success : bool
855 True if the algorithm converged to the desired accuracy (``status=0``).
857 Notes
858 -----
859 This function implements a 4th order collocation algorithm with the
860 control of residuals similar to [1]_. A collocation system is solved
861 by a damped Newton method with an affine-invariant criterion function as
862 described in [3]_.
864 Note that in [1]_ integral residuals are defined without normalization
865 by interval lengths. So, their definition is different by a multiplier of
866 h**0.5 (h is an interval length) from the definition used here.
868 .. versionadded:: 0.18.0
870 References
871 ----------
872 .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
873 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
874 Number 3, pp. 299-316, 2001.
875 .. [2] L.F. Shampine, P. H. Muir and H. Xu, "A User-Friendly Fortran BVP
876 Solver".
877 .. [3] U. Ascher, R. Mattheij and R. Russell "Numerical Solution of
878 Boundary Value Problems for Ordinary Differential Equations".
879 .. [4] `Cauchy-Riemann equations
880 <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
881 Wikipedia.
883 Examples
884 --------
885 In the first example, we solve Bratu's problem::
887 y'' + k * exp(y) = 0
888 y(0) = y(1) = 0
890 for k = 1.
892 We rewrite the equation as a first-order system and implement its
893 right-hand side evaluation::
895 y1' = y2
896 y2' = -exp(y1)
898 >>> import numpy as np
899 >>> def fun(x, y):
900 ... return np.vstack((y[1], -np.exp(y[0])))
902 Implement evaluation of the boundary condition residuals:
904 >>> def bc(ya, yb):
905 ... return np.array([ya[0], yb[0]])
907 Define the initial mesh with 5 nodes:
909 >>> x = np.linspace(0, 1, 5)
911 This problem is known to have two solutions. To obtain both of them, we
912 use two different initial guesses for y. We denote them by subscripts
913 a and b.
915 >>> y_a = np.zeros((2, x.size))
916 >>> y_b = np.zeros((2, x.size))
917 >>> y_b[0] = 3
919 Now we are ready to run the solver.
921 >>> from scipy.integrate import solve_bvp
922 >>> res_a = solve_bvp(fun, bc, x, y_a)
923 >>> res_b = solve_bvp(fun, bc, x, y_b)
925 Let's plot the two found solutions. We take an advantage of having the
926 solution in a spline form to produce a smooth plot.
928 >>> x_plot = np.linspace(0, 1, 100)
929 >>> y_plot_a = res_a.sol(x_plot)[0]
930 >>> y_plot_b = res_b.sol(x_plot)[0]
931 >>> import matplotlib.pyplot as plt
932 >>> plt.plot(x_plot, y_plot_a, label='y_a')
933 >>> plt.plot(x_plot, y_plot_b, label='y_b')
934 >>> plt.legend()
935 >>> plt.xlabel("x")
936 >>> plt.ylabel("y")
937 >>> plt.show()
939 We see that the two solutions have similar shape, but differ in scale
940 significantly.
942 In the second example, we solve a simple Sturm-Liouville problem::
944 y'' + k**2 * y = 0
945 y(0) = y(1) = 0
947 It is known that a non-trivial solution y = A * sin(k * x) is possible for
948 k = pi * n, where n is an integer. To establish the normalization constant
949 A = 1 we add a boundary condition::
951 y'(0) = k
953 Again, we rewrite our equation as a first-order system and implement its
954 right-hand side evaluation::
956 y1' = y2
957 y2' = -k**2 * y1
959 >>> def fun(x, y, p):
960 ... k = p[0]
961 ... return np.vstack((y[1], -k**2 * y[0]))
963 Note that parameters p are passed as a vector (with one element in our
964 case).
966 Implement the boundary conditions:
968 >>> def bc(ya, yb, p):
969 ... k = p[0]
970 ... return np.array([ya[0], yb[0], ya[1] - k])
972 Set up the initial mesh and guess for y. We aim to find the solution for
973 k = 2 * pi, to achieve that we set values of y to approximately follow
974 sin(2 * pi * x):
976 >>> x = np.linspace(0, 1, 5)
977 >>> y = np.zeros((2, x.size))
978 >>> y[0, 1] = 1
979 >>> y[0, 3] = -1
981 Run the solver with 6 as an initial guess for k.
983 >>> sol = solve_bvp(fun, bc, x, y, p=[6])
985 We see that the found k is approximately correct:
987 >>> sol.p[0]
988 6.28329460046
990 And, finally, plot the solution to see the anticipated sinusoid:
992 >>> x_plot = np.linspace(0, 1, 100)
993 >>> y_plot = sol.sol(x_plot)[0]
994 >>> plt.plot(x_plot, y_plot)
995 >>> plt.xlabel("x")
996 >>> plt.ylabel("y")
997 >>> plt.show()
998 """
999 x = np.asarray(x, dtype=float)
1000 if x.ndim != 1:
1001 raise ValueError("`x` must be 1 dimensional.")
1002 h = np.diff(x)
1003 if np.any(h <= 0):
1004 raise ValueError("`x` must be strictly increasing.")
1005 a = x[0]
1007 y = np.asarray(y)
1008 if np.issubdtype(y.dtype, np.complexfloating):
1009 dtype = complex
1010 else:
1011 dtype = float
1012 y = y.astype(dtype, copy=False)
1014 if y.ndim != 2:
1015 raise ValueError("`y` must be 2 dimensional.")
1016 if y.shape[1] != x.shape[0]:
1017 raise ValueError("`y` is expected to have {} columns, but actually "
1018 "has {}.".format(x.shape[0], y.shape[1]))
1020 if p is None:
1021 p = np.array([])
1022 else:
1023 p = np.asarray(p, dtype=dtype)
1024 if p.ndim != 1:
1025 raise ValueError("`p` must be 1 dimensional.")
1027 if tol < 100 * EPS:
1028 warn("`tol` is too low, setting to {:.2e}".format(100 * EPS))
1029 tol = 100 * EPS
1031 if verbose not in [0, 1, 2]:
1032 raise ValueError("`verbose` must be in [0, 1, 2].")
1034 n = y.shape[0]
1035 k = p.shape[0]
1037 if S is not None:
1038 S = np.asarray(S, dtype=dtype)
1039 if S.shape != (n, n):
1040 raise ValueError("`S` is expected to have shape {}, "
1041 "but actually has {}".format((n, n), S.shape))
1043 # Compute I - S^+ S to impose necessary boundary conditions.
1044 B = np.identity(n) - np.dot(pinv(S), S)
1046 y[:, 0] = np.dot(B, y[:, 0])
1048 # Compute (I - S)^+ to correct derivatives at x=a.
1049 D = pinv(np.identity(n) - S)
1050 else:
1051 B = None
1052 D = None
1054 if bc_tol is None:
1055 bc_tol = tol
1057 # Maximum number of iterations
1058 max_iteration = 10
1060 fun_wrapped, bc_wrapped, fun_jac_wrapped, bc_jac_wrapped = wrap_functions(
1061 fun, bc, fun_jac, bc_jac, k, a, S, D, dtype)
1063 f = fun_wrapped(x, y, p)
1064 if f.shape != y.shape:
1065 raise ValueError("`fun` return is expected to have shape {}, "
1066 "but actually has {}.".format(y.shape, f.shape))
1068 bc_res = bc_wrapped(y[:, 0], y[:, -1], p)
1069 if bc_res.shape != (n + k,):
1070 raise ValueError("`bc` return is expected to have shape {}, "
1071 "but actually has {}.".format((n + k,), bc_res.shape))
1073 status = 0
1074 iteration = 0
1075 if verbose == 2:
1076 print_iteration_header()
1078 while True:
1079 m = x.shape[0]
1081 col_fun, jac_sys = prepare_sys(n, m, k, fun_wrapped, bc_wrapped,
1082 fun_jac_wrapped, bc_jac_wrapped, x, h)
1083 y, p, singular = solve_newton(n, m, h, col_fun, bc_wrapped, jac_sys,
1084 y, p, B, tol, bc_tol)
1085 iteration += 1
1087 col_res, y_middle, f, f_middle = collocation_fun(fun_wrapped, y,
1088 p, x, h)
1089 bc_res = bc_wrapped(y[:, 0], y[:, -1], p)
1090 max_bc_res = np.max(abs(bc_res))
1092 # This relation is not trivial, but can be verified.
1093 r_middle = 1.5 * col_res / h
1094 sol = create_spline(y, f, x, h)
1095 rms_res = estimate_rms_residuals(fun_wrapped, sol, x, h, p,
1096 r_middle, f_middle)
1097 max_rms_res = np.max(rms_res)
1099 if singular:
1100 status = 2
1101 break
1103 insert_1, = np.nonzero((rms_res > tol) & (rms_res < 100 * tol))
1104 insert_2, = np.nonzero(rms_res >= 100 * tol)
1105 nodes_added = insert_1.shape[0] + 2 * insert_2.shape[0]
1107 if m + nodes_added > max_nodes:
1108 status = 1
1109 if verbose == 2:
1110 nodes_added = "({})".format(nodes_added)
1111 print_iteration_progress(iteration, max_rms_res, max_bc_res,
1112 m, nodes_added)
1113 break
1115 if verbose == 2:
1116 print_iteration_progress(iteration, max_rms_res, max_bc_res, m,
1117 nodes_added)
1119 if nodes_added > 0:
1120 x = modify_mesh(x, insert_1, insert_2)
1121 h = np.diff(x)
1122 y = sol(x)
1123 elif max_bc_res <= bc_tol:
1124 status = 0
1125 break
1126 elif iteration >= max_iteration:
1127 status = 3
1128 break
1130 if verbose > 0:
1131 if status == 0:
1132 print("Solved in {} iterations, number of nodes {}. \n"
1133 "Maximum relative residual: {:.2e} \n"
1134 "Maximum boundary residual: {:.2e}"
1135 .format(iteration, x.shape[0], max_rms_res, max_bc_res))
1136 elif status == 1:
1137 print("Number of nodes is exceeded after iteration {}. \n"
1138 "Maximum relative residual: {:.2e} \n"
1139 "Maximum boundary residual: {:.2e}"
1140 .format(iteration, max_rms_res, max_bc_res))
1141 elif status == 2:
1142 print("Singular Jacobian encountered when solving the collocation "
1143 "system on iteration {}. \n"
1144 "Maximum relative residual: {:.2e} \n"
1145 "Maximum boundary residual: {:.2e}"
1146 .format(iteration, max_rms_res, max_bc_res))
1147 elif status == 3:
1148 print("The solver was unable to satisfy boundary conditions "
1149 "tolerance on iteration {}. \n"
1150 "Maximum relative residual: {:.2e} \n"
1151 "Maximum boundary residual: {:.2e}"
1152 .format(iteration, max_rms_res, max_bc_res))
1154 if p.size == 0:
1155 p = None
1157 return BVPResult(sol=sol, p=p, x=x, y=y, yp=f, rms_residuals=rms_res,
1158 niter=iteration, status=status,
1159 message=TERMINATION_MESSAGES[status], success=status == 0)