Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/integrate/_bvp.py : 7%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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 In our case, a[i, :, :] and b[i, :, :] are always square.
147 """
148 # Empirical optimization. Use outer Python loop and BLAS for large
149 # matrices, otherwise use a single einsum call.
150 if a.shape[1] > 50:
151 out = np.empty_like(a)
152 for i in range(a.shape[0]):
153 out[i] = np.dot(a[i], b[i])
154 return out
155 else:
156 return np.einsum('...ij,...jk->...ik', a, b)
159def construct_global_jac(n, m, k, i_jac, j_jac, h, df_dy, df_dy_middle, df_dp,
160 df_dp_middle, dbc_dya, dbc_dyb, dbc_dp):
161 """Construct the Jacobian of the collocation system.
163 There are n * m + k functions: m - 1 collocations residuals, each
164 containing n components, followed by n + k boundary condition residuals.
166 There are n * m + k variables: m vectors of y, each containing n
167 components, followed by k values of vector p.
169 For example, let m = 4, n = 2 and k = 1, then the Jacobian will have
170 the following sparsity structure:
172 1 1 2 2 0 0 0 0 5
173 1 1 2 2 0 0 0 0 5
174 0 0 1 1 2 2 0 0 5
175 0 0 1 1 2 2 0 0 5
176 0 0 0 0 1 1 2 2 5
177 0 0 0 0 1 1 2 2 5
179 3 3 0 0 0 0 4 4 6
180 3 3 0 0 0 0 4 4 6
181 3 3 0 0 0 0 4 4 6
183 Zeros denote identically zero values, other values denote different kinds
184 of blocks in the matrix (see below). The blank row indicates the separation
185 of collocation residuals from boundary conditions. And the blank column
186 indicates the separation of y values from p values.
188 Refer to [1]_ (p. 306) for the formula of n x n blocks for derivatives
189 of collocation residuals with respect to y.
191 Parameters
192 ----------
193 n : int
194 Number of equations in the ODE system.
195 m : int
196 Number of nodes in the mesh.
197 k : int
198 Number of the unknown parameters.
199 i_jac, j_jac : ndarray
200 Row and column indices returned by `compute_jac_indices`. They
201 represent different blocks in the Jacobian matrix in the following
202 order (see the scheme above):
204 * 1: m - 1 diagonal n x n blocks for the collocation residuals.
205 * 2: m - 1 off-diagonal n x n blocks for the collocation residuals.
206 * 3 : (n + k) x n block for the dependency of the boundary
207 conditions on ya.
208 * 4: (n + k) x n block for the dependency of the boundary
209 conditions on yb.
210 * 5: (m - 1) * n x k block for the dependency of the collocation
211 residuals on p.
212 * 6: (n + k) x k block for the dependency of the boundary
213 conditions on p.
215 df_dy : ndarray, shape (n, n, m)
216 Jacobian of f with respect to y computed at the mesh nodes.
217 df_dy_middle : ndarray, shape (n, n, m - 1)
218 Jacobian of f with respect to y computed at the middle between the
219 mesh nodes.
220 df_dp : ndarray with shape (n, k, m) or None
221 Jacobian of f with respect to p computed at the mesh nodes.
222 df_dp_middle: ndarray with shape (n, k, m - 1) or None
223 Jacobian of f with respect to p computed at the middle between the
224 mesh nodes.
225 dbc_dya, dbc_dyb : ndarray, shape (n, n)
226 Jacobian of bc with respect to ya and yb.
227 dbc_dp: ndarray with shape (n, k) or None
228 Jacobian of bc with respect to p.
230 Returns
231 -------
232 J : csc_matrix, shape (n * m + k, n * m + k)
233 Jacobian of the collocation system in a sparse form.
235 References
236 ----------
237 .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
238 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
239 Number 3, pp. 299-316, 2001.
240 """
241 df_dy = np.transpose(df_dy, (2, 0, 1))
242 df_dy_middle = np.transpose(df_dy_middle, (2, 0, 1))
244 h = h[:, np.newaxis, np.newaxis]
246 dtype = df_dy.dtype
248 # Computing diagonal n x n blocks.
249 dPhi_dy_0 = np.empty((m - 1, n, n), dtype=dtype)
250 dPhi_dy_0[:] = -np.identity(n)
251 dPhi_dy_0 -= h / 6 * (df_dy[:-1] + 2 * df_dy_middle)
252 T = stacked_matmul(df_dy_middle, df_dy[:-1])
253 dPhi_dy_0 -= h**2 / 12 * T
255 # Computing off-diagonal n x n blocks.
256 dPhi_dy_1 = np.empty((m - 1, n, n), dtype=dtype)
257 dPhi_dy_1[:] = np.identity(n)
258 dPhi_dy_1 -= h / 6 * (df_dy[1:] + 2 * df_dy_middle)
259 T = stacked_matmul(df_dy_middle, df_dy[1:])
260 dPhi_dy_1 += h**2 / 12 * T
262 values = np.hstack((dPhi_dy_0.ravel(), dPhi_dy_1.ravel(), dbc_dya.ravel(),
263 dbc_dyb.ravel()))
265 if k > 0:
266 df_dp = np.transpose(df_dp, (2, 0, 1))
267 df_dp_middle = np.transpose(df_dp_middle, (2, 0, 1))
268 T = stacked_matmul(df_dy_middle, df_dp[:-1] - df_dp[1:])
269 df_dp_middle += 0.125 * h * T
270 dPhi_dp = -h/6 * (df_dp[:-1] + df_dp[1:] + 4 * df_dp_middle)
271 values = np.hstack((values, dPhi_dp.ravel(), dbc_dp.ravel()))
273 J = coo_matrix((values, (i_jac, j_jac)))
274 return csc_matrix(J)
277def collocation_fun(fun, y, p, x, h):
278 """Evaluate collocation residuals.
280 This function lies in the core of the method. The solution is sought
281 as a cubic C1 continuous spline with derivatives matching the ODE rhs
282 at given nodes `x`. Collocation conditions are formed from the equality
283 of the spline derivatives and rhs of the ODE system in the middle points
284 between nodes.
286 Such method is classified to Lobbato IIIA family in ODE literature.
287 Refer to [1]_ for the formula and some discussion.
289 Returns
290 -------
291 col_res : ndarray, shape (n, m - 1)
292 Collocation residuals at the middle points of the mesh intervals.
293 y_middle : ndarray, shape (n, m - 1)
294 Values of the cubic spline evaluated at the middle points of the mesh
295 intervals.
296 f : ndarray, shape (n, m)
297 RHS of the ODE system evaluated at the mesh nodes.
298 f_middle : ndarray, shape (n, m - 1)
299 RHS of the ODE system evaluated at the middle points of the mesh
300 intervals (and using `y_middle`).
302 References
303 ----------
304 .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
305 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
306 Number 3, pp. 299-316, 2001.
307 """
308 f = fun(x, y, p)
309 y_middle = (0.5 * (y[:, 1:] + y[:, :-1]) -
310 0.125 * h * (f[:, 1:] - f[:, :-1]))
311 f_middle = fun(x[:-1] + 0.5 * h, y_middle, p)
312 col_res = y[:, 1:] - y[:, :-1] - h / 6 * (f[:, :-1] + f[:, 1:] +
313 4 * f_middle)
315 return col_res, y_middle, f, f_middle
318def prepare_sys(n, m, k, fun, bc, fun_jac, bc_jac, x, h):
319 """Create the function and the Jacobian for the collocation system."""
320 x_middle = x[:-1] + 0.5 * h
321 i_jac, j_jac = compute_jac_indices(n, m, k)
323 def col_fun(y, p):
324 return collocation_fun(fun, y, p, x, h)
326 def sys_jac(y, p, y_middle, f, f_middle, bc0):
327 if fun_jac is None:
328 df_dy, df_dp = estimate_fun_jac(fun, x, y, p, f)
329 df_dy_middle, df_dp_middle = estimate_fun_jac(
330 fun, x_middle, y_middle, p, f_middle)
331 else:
332 df_dy, df_dp = fun_jac(x, y, p)
333 df_dy_middle, df_dp_middle = fun_jac(x_middle, y_middle, p)
335 if bc_jac is None:
336 dbc_dya, dbc_dyb, dbc_dp = estimate_bc_jac(bc, y[:, 0], y[:, -1],
337 p, bc0)
338 else:
339 dbc_dya, dbc_dyb, dbc_dp = bc_jac(y[:, 0], y[:, -1], p)
341 return construct_global_jac(n, m, k, i_jac, j_jac, h, df_dy,
342 df_dy_middle, df_dp, df_dp_middle, dbc_dya,
343 dbc_dyb, dbc_dp)
345 return col_fun, sys_jac
348def solve_newton(n, m, h, col_fun, bc, jac, y, p, B, bvp_tol, bc_tol):
349 """Solve the nonlinear collocation system by a Newton method.
351 This is a simple Newton method with a backtracking line search. As
352 advised in [1]_, an affine-invariant criterion function F = ||J^-1 r||^2
353 is used, where J is the Jacobian matrix at the current iteration and r is
354 the vector or collocation residuals (values of the system lhs).
356 The method alters between full Newton iterations and the fixed-Jacobian
357 iterations based
359 There are other tricks proposed in [1]_, but they are not used as they
360 don't seem to improve anything significantly, and even break the
361 convergence on some test problems I tried.
363 All important parameters of the algorithm are defined inside the function.
365 Parameters
366 ----------
367 n : int
368 Number of equations in the ODE system.
369 m : int
370 Number of nodes in the mesh.
371 h : ndarray, shape (m-1,)
372 Mesh intervals.
373 col_fun : callable
374 Function computing collocation residuals.
375 bc : callable
376 Function computing boundary condition residuals.
377 jac : callable
378 Function computing the Jacobian of the whole system (including
379 collocation and boundary condition residuals). It is supposed to
380 return csc_matrix.
381 y : ndarray, shape (n, m)
382 Initial guess for the function values at the mesh nodes.
383 p : ndarray, shape (k,)
384 Initial guess for the unknown parameters.
385 B : ndarray with shape (n, n) or None
386 Matrix to force the S y(a) = 0 condition for a problems with the
387 singular term. If None, the singular term is assumed to be absent.
388 bvp_tol : float
389 Tolerance to which we want to solve a BVP.
390 bc_tol : float
391 Tolerance to which we want to satisfy the boundary conditions.
393 Returns
394 -------
395 y : ndarray, shape (n, m)
396 Final iterate for the function values at the mesh nodes.
397 p : ndarray, shape (k,)
398 Final iterate for the unknown parameters.
399 singular : bool
400 True, if the LU decomposition failed because Jacobian turned out
401 to be singular.
403 References
404 ----------
405 .. [1] U. Ascher, R. Mattheij and R. Russell "Numerical Solution of
406 Boundary Value Problems for Ordinary Differential Equations"
407 """
408 # We know that the solution residuals at the middle points of the mesh
409 # are connected with collocation residuals r_middle = 1.5 * col_res / h.
410 # As our BVP solver tries to decrease relative residuals below a certain
411 # tolerance, it seems reasonable to terminated Newton iterations by
412 # comparison of r_middle / (1 + np.abs(f_middle)) with a certain threshold,
413 # which we choose to be 1.5 orders lower than the BVP tolerance. We rewrite
414 # the condition as col_res < tol_r * (1 + np.abs(f_middle)), then tol_r
415 # should be computed as follows:
416 tol_r = 2/3 * h * 5e-2 * bvp_tol
418 # Maximum allowed number of Jacobian evaluation and factorization, in
419 # other words, the maximum number of full Newton iterations. A small value
420 # is recommended in the literature.
421 max_njev = 4
423 # Maximum number of iterations, considering that some of them can be
424 # performed with the fixed Jacobian. In theory, such iterations are cheap,
425 # but it's not that simple in Python.
426 max_iter = 8
428 # Minimum relative improvement of the criterion function to accept the
429 # step (Armijo constant).
430 sigma = 0.2
432 # Step size decrease factor for backtracking.
433 tau = 0.5
435 # Maximum number of backtracking steps, the minimum step is then
436 # tau ** n_trial.
437 n_trial = 4
439 col_res, y_middle, f, f_middle = col_fun(y, p)
440 bc_res = bc(y[:, 0], y[:, -1], p)
441 res = np.hstack((col_res.ravel(order='F'), bc_res))
443 njev = 0
444 singular = False
445 recompute_jac = True
446 for iteration in range(max_iter):
447 if recompute_jac:
448 J = jac(y, p, y_middle, f, f_middle, bc_res)
449 njev += 1
450 try:
451 LU = splu(J)
452 except RuntimeError:
453 singular = True
454 break
456 step = LU.solve(res)
457 cost = np.dot(step, step)
459 y_step = step[:m * n].reshape((n, m), order='F')
460 p_step = step[m * n:]
462 alpha = 1
463 for trial in range(n_trial + 1):
464 y_new = y - alpha * y_step
465 if B is not None:
466 y_new[:, 0] = np.dot(B, y_new[:, 0])
467 p_new = p - alpha * p_step
469 col_res, y_middle, f, f_middle = col_fun(y_new, p_new)
470 bc_res = bc(y_new[:, 0], y_new[:, -1], p_new)
471 res = np.hstack((col_res.ravel(order='F'), bc_res))
473 step_new = LU.solve(res)
474 cost_new = np.dot(step_new, step_new)
475 if cost_new < (1 - 2 * alpha * sigma) * cost:
476 break
478 if trial < n_trial:
479 alpha *= tau
481 y = y_new
482 p = p_new
484 if njev == max_njev:
485 break
487 if (np.all(np.abs(col_res) < tol_r * (1 + np.abs(f_middle))) and
488 np.all(np.abs(bc_res) < bc_tol)):
489 break
491 # If the full step was taken, then we are going to continue with
492 # the same Jacobian. This is the approach of BVP_SOLVER.
493 if alpha == 1:
494 step = step_new
495 cost = cost_new
496 recompute_jac = False
497 else:
498 recompute_jac = True
500 return y, p, singular
503def print_iteration_header():
504 print("{:^15}{:^15}{:^15}{:^15}{:^15}".format(
505 "Iteration", "Max residual", "Max BC residual", "Total nodes",
506 "Nodes added"))
509def print_iteration_progress(iteration, residual, bc_residual, total_nodes,
510 nodes_added):
511 print("{:^15}{:^15.2e}{:^15.2e}{:^15}{:^15}".format(
512 iteration, residual, bc_residual, total_nodes, nodes_added))
515class BVPResult(OptimizeResult):
516 pass
519TERMINATION_MESSAGES = {
520 0: "The algorithm converged to the desired accuracy.",
521 1: "The maximum number of mesh nodes is exceeded.",
522 2: "A singular Jacobian encountered when solving the collocation system.",
523 3: "The solver was unable to satisfy boundary conditions tolerance on iteration 10."
524}
527def estimate_rms_residuals(fun, sol, x, h, p, r_middle, f_middle):
528 """Estimate rms values of collocation residuals using Lobatto quadrature.
530 The residuals are defined as the difference between the derivatives of
531 our solution and rhs of the ODE system. We use relative residuals, i.e.,
532 normalized by 1 + np.abs(f). RMS values are computed as sqrt from the
533 normalized integrals of the squared relative residuals over each interval.
534 Integrals are estimated using 5-point Lobatto quadrature [1]_, we use the
535 fact that residuals at the mesh nodes are identically zero.
537 In [2] they don't normalize integrals by interval lengths, which gives
538 a higher rate of convergence of the residuals by the factor of h**0.5.
539 I chose to do such normalization for an ease of interpretation of return
540 values as RMS estimates.
542 Returns
543 -------
544 rms_res : ndarray, shape (m - 1,)
545 Estimated rms values of the relative residuals over each interval.
547 References
548 ----------
549 .. [1] http://mathworld.wolfram.com/LobattoQuadrature.html
550 .. [2] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
551 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
552 Number 3, pp. 299-316, 2001.
553 """
554 x_middle = x[:-1] + 0.5 * h
555 s = 0.5 * h * (3/7)**0.5
556 x1 = x_middle + s
557 x2 = x_middle - s
558 y1 = sol(x1)
559 y2 = sol(x2)
560 y1_prime = sol(x1, 1)
561 y2_prime = sol(x2, 1)
562 f1 = fun(x1, y1, p)
563 f2 = fun(x2, y2, p)
564 r1 = y1_prime - f1
565 r2 = y2_prime - f2
567 r_middle /= 1 + np.abs(f_middle)
568 r1 /= 1 + np.abs(f1)
569 r2 /= 1 + np.abs(f2)
571 r1 = np.sum(np.real(r1 * np.conj(r1)), axis=0)
572 r2 = np.sum(np.real(r2 * np.conj(r2)), axis=0)
573 r_middle = np.sum(np.real(r_middle * np.conj(r_middle)), axis=0)
575 return (0.5 * (32 / 45 * r_middle + 49 / 90 * (r1 + r2))) ** 0.5
578def create_spline(y, yp, x, h):
579 """Create a cubic spline given values and derivatives.
581 Formulas for the coefficients are taken from interpolate.CubicSpline.
583 Returns
584 -------
585 sol : PPoly
586 Constructed spline as a PPoly instance.
587 """
588 from scipy.interpolate import PPoly
590 n, m = y.shape
591 c = np.empty((4, n, m - 1), dtype=y.dtype)
592 slope = (y[:, 1:] - y[:, :-1]) / h
593 t = (yp[:, :-1] + yp[:, 1:] - 2 * slope) / h
594 c[0] = t / h
595 c[1] = (slope - yp[:, :-1]) / h - t
596 c[2] = yp[:, :-1]
597 c[3] = y[:, :-1]
598 c = np.rollaxis(c, 1)
600 return PPoly(c, x, extrapolate=True, axis=1)
603def modify_mesh(x, insert_1, insert_2):
604 """Insert nodes into a mesh.
606 Nodes removal logic is not established, its impact on the solver is
607 presumably negligible. So, only insertion is done in this function.
609 Parameters
610 ----------
611 x : ndarray, shape (m,)
612 Mesh nodes.
613 insert_1 : ndarray
614 Intervals to each insert 1 new node in the middle.
615 insert_2 : ndarray
616 Intervals to each insert 2 new nodes, such that divide an interval
617 into 3 equal parts.
619 Returns
620 -------
621 x_new : ndarray
622 New mesh nodes.
624 Notes
625 -----
626 `insert_1` and `insert_2` should not have common values.
627 """
628 # Because np.insert implementation apparently varies with a version of
629 # NumPy, we use a simple and reliable approach with sorting.
630 return np.sort(np.hstack((
631 x,
632 0.5 * (x[insert_1] + x[insert_1 + 1]),
633 (2 * x[insert_2] + x[insert_2 + 1]) / 3,
634 (x[insert_2] + 2 * x[insert_2 + 1]) / 3
635 )))
638def wrap_functions(fun, bc, fun_jac, bc_jac, k, a, S, D, dtype):
639 """Wrap functions for unified usage in the solver."""
640 if fun_jac is None:
641 fun_jac_wrapped = None
643 if bc_jac is None:
644 bc_jac_wrapped = None
646 if k == 0:
647 def fun_p(x, y, _):
648 return np.asarray(fun(x, y), dtype)
650 def bc_wrapped(ya, yb, _):
651 return np.asarray(bc(ya, yb), dtype)
653 if fun_jac is not None:
654 def fun_jac_p(x, y, _):
655 return np.asarray(fun_jac(x, y), dtype), None
657 if bc_jac is not None:
658 def bc_jac_wrapped(ya, yb, _):
659 dbc_dya, dbc_dyb = bc_jac(ya, yb)
660 return (np.asarray(dbc_dya, dtype),
661 np.asarray(dbc_dyb, dtype), None)
662 else:
663 def fun_p(x, y, p):
664 return np.asarray(fun(x, y, p), dtype)
666 def bc_wrapped(x, y, p):
667 return np.asarray(bc(x, y, p), dtype)
669 if fun_jac is not None:
670 def fun_jac_p(x, y, p):
671 df_dy, df_dp = fun_jac(x, y, p)
672 return np.asarray(df_dy, dtype), np.asarray(df_dp, dtype)
674 if bc_jac is not None:
675 def bc_jac_wrapped(ya, yb, p):
676 dbc_dya, dbc_dyb, dbc_dp = bc_jac(ya, yb, p)
677 return (np.asarray(dbc_dya, dtype), np.asarray(dbc_dyb, dtype),
678 np.asarray(dbc_dp, dtype))
680 if S is None:
681 fun_wrapped = fun_p
682 else:
683 def fun_wrapped(x, y, p):
684 f = fun_p(x, y, p)
685 if x[0] == a:
686 f[:, 0] = np.dot(D, f[:, 0])
687 f[:, 1:] += np.dot(S, y[:, 1:]) / (x[1:] - a)
688 else:
689 f += np.dot(S, y) / (x - a)
690 return f
692 if fun_jac is not None:
693 if S is None:
694 fun_jac_wrapped = fun_jac_p
695 else:
696 Sr = S[:, :, np.newaxis]
698 def fun_jac_wrapped(x, y, p):
699 df_dy, df_dp = fun_jac_p(x, y, p)
700 if x[0] == a:
701 df_dy[:, :, 0] = np.dot(D, df_dy[:, :, 0])
702 df_dy[:, :, 1:] += Sr / (x[1:] - a)
703 else:
704 df_dy += Sr / (x - a)
706 return df_dy, df_dp
708 return fun_wrapped, bc_wrapped, fun_jac_wrapped, bc_jac_wrapped
711def solve_bvp(fun, bc, x, y, p=None, S=None, fun_jac=None, bc_jac=None,
712 tol=1e-3, max_nodes=1000, verbose=0, bc_tol=None):
713 """Solve a boundary value problem for a system of ODEs.
715 This function numerically solves a first order system of ODEs subject to
716 two-point boundary conditions::
718 dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b
719 bc(y(a), y(b), p) = 0
721 Here x is a 1-D independent variable, y(x) is an N-D
722 vector-valued function and p is a k-D vector of unknown
723 parameters which is to be found along with y(x). For the problem to be
724 determined, there must be n + k boundary conditions, i.e., bc must be an
725 (n + k)-D function.
727 The last singular term on the right-hand side of the system is optional.
728 It is defined by an n-by-n matrix S, such that the solution must satisfy
729 S y(a) = 0. This condition will be forced during iterations, so it must not
730 contradict boundary conditions. See [2]_ for the explanation how this term
731 is handled when solving BVPs numerically.
733 Problems in a complex domain can be solved as well. In this case, y and p
734 are considered to be complex, and f and bc are assumed to be complex-valued
735 functions, but x stays real. Note that f and bc must be complex
736 differentiable (satisfy Cauchy-Riemann equations [4]_), otherwise you
737 should rewrite your problem for real and imaginary parts separately. To
738 solve a problem in a complex domain, pass an initial guess for y with a
739 complex data type (see below).
741 Parameters
742 ----------
743 fun : callable
744 Right-hand side of the system. The calling signature is ``fun(x, y)``,
745 or ``fun(x, y, p)`` if parameters are present. All arguments are
746 ndarray: ``x`` with shape (m,), ``y`` with shape (n, m), meaning that
747 ``y[:, i]`` corresponds to ``x[i]``, and ``p`` with shape (k,). The
748 return value must be an array with shape (n, m) and with the same
749 layout as ``y``.
750 bc : callable
751 Function evaluating residuals of the boundary conditions. The calling
752 signature is ``bc(ya, yb)``, or ``bc(ya, yb, p)`` if parameters are
753 present. All arguments are ndarray: ``ya`` and ``yb`` with shape (n,),
754 and ``p`` with shape (k,). The return value must be an array with
755 shape (n + k,).
756 x : array_like, shape (m,)
757 Initial mesh. Must be a strictly increasing sequence of real numbers
758 with ``x[0]=a`` and ``x[-1]=b``.
759 y : array_like, shape (n, m)
760 Initial guess for the function values at the mesh nodes, ith column
761 corresponds to ``x[i]``. For problems in a complex domain pass `y`
762 with a complex data type (even if the initial guess is purely real).
763 p : array_like with shape (k,) or None, optional
764 Initial guess for the unknown parameters. If None (default), it is
765 assumed that the problem doesn't depend on any parameters.
766 S : array_like with shape (n, n) or None
767 Matrix defining the singular term. If None (default), the problem is
768 solved without the singular term.
769 fun_jac : callable or None, optional
770 Function computing derivatives of f with respect to y and p. The
771 calling signature is ``fun_jac(x, y)``, or ``fun_jac(x, y, p)`` if
772 parameters are present. The return must contain 1 or 2 elements in the
773 following order:
775 * df_dy : array_like with shape (n, n, m), where an element
776 (i, j, q) equals to d f_i(x_q, y_q, p) / d (y_q)_j.
777 * df_dp : array_like with shape (n, k, m), where an element
778 (i, j, q) equals to d f_i(x_q, y_q, p) / d p_j.
780 Here q numbers nodes at which x and y are defined, whereas i and j
781 number vector components. If the problem is solved without unknown
782 parameters, df_dp should not be returned.
784 If `fun_jac` is None (default), the derivatives will be estimated
785 by the forward finite differences.
786 bc_jac : callable or None, optional
787 Function computing derivatives of bc with respect to ya, yb, and p.
788 The calling signature is ``bc_jac(ya, yb)``, or ``bc_jac(ya, yb, p)``
789 if parameters are present. The return must contain 2 or 3 elements in
790 the following order:
792 * dbc_dya : array_like with shape (n, n), where an element (i, j)
793 equals to d bc_i(ya, yb, p) / d ya_j.
794 * dbc_dyb : array_like with shape (n, n), where an element (i, j)
795 equals to d bc_i(ya, yb, p) / d yb_j.
796 * dbc_dp : array_like with shape (n, k), where an element (i, j)
797 equals to d bc_i(ya, yb, p) / d p_j.
799 If the problem is solved without unknown parameters, dbc_dp should not
800 be returned.
802 If `bc_jac` is None (default), the derivatives will be estimated by
803 the forward finite differences.
804 tol : float, optional
805 Desired tolerance of the solution. If we define ``r = y' - f(x, y)``,
806 where y is the found solution, then the solver tries to achieve on each
807 mesh interval ``norm(r / (1 + abs(f)) < tol``, where ``norm`` is
808 estimated in a root mean squared sense (using a numerical quadrature
809 formula). Default is 1e-3.
810 max_nodes : int, optional
811 Maximum allowed number of the mesh nodes. If exceeded, the algorithm
812 terminates. Default is 1000.
813 verbose : {0, 1, 2}, optional
814 Level of algorithm's verbosity:
816 * 0 (default) : work silently.
817 * 1 : display a termination report.
818 * 2 : display progress during iterations.
819 bc_tol : float, optional
820 Desired absolute tolerance for the boundary condition residuals: `bc`
821 value should satisfy ``abs(bc) < bc_tol`` component-wise.
822 Equals to `tol` by default. Up to 10 iterations are allowed to achieve this
823 tolerance.
825 Returns
826 -------
827 Bunch object with the following fields defined:
828 sol : PPoly
829 Found solution for y as `scipy.interpolate.PPoly` instance, a C1
830 continuous cubic spline.
831 p : ndarray or None, shape (k,)
832 Found parameters. None, if the parameters were not present in the
833 problem.
834 x : ndarray, shape (m,)
835 Nodes of the final mesh.
836 y : ndarray, shape (n, m)
837 Solution values at the mesh nodes.
838 yp : ndarray, shape (n, m)
839 Solution derivatives at the mesh nodes.
840 rms_residuals : ndarray, shape (m - 1,)
841 RMS values of the relative residuals over each mesh interval (see the
842 description of `tol` parameter).
843 niter : int
844 Number of completed iterations.
845 status : int
846 Reason for algorithm termination:
848 * 0: The algorithm converged to the desired accuracy.
849 * 1: The maximum number of mesh nodes is exceeded.
850 * 2: A singular Jacobian encountered when solving the collocation
851 system.
853 message : string
854 Verbal description of the termination reason.
855 success : bool
856 True if the algorithm converged to the desired accuracy (``status=0``).
858 Notes
859 -----
860 This function implements a 4th order collocation algorithm with the
861 control of residuals similar to [1]_. A collocation system is solved
862 by a damped Newton method with an affine-invariant criterion function as
863 described in [3]_.
865 Note that in [1]_ integral residuals are defined without normalization
866 by interval lengths. So, their definition is different by a multiplier of
867 h**0.5 (h is an interval length) from the definition used here.
869 .. versionadded:: 0.18.0
871 References
872 ----------
873 .. [1] J. Kierzenka, L. F. Shampine, "A BVP Solver Based on Residual
874 Control and the Maltab PSE", ACM Trans. Math. Softw., Vol. 27,
875 Number 3, pp. 299-316, 2001.
876 .. [2] L.F. Shampine, P. H. Muir and H. Xu, "A User-Friendly Fortran BVP
877 Solver".
878 .. [3] U. Ascher, R. Mattheij and R. Russell "Numerical Solution of
879 Boundary Value Problems for Ordinary Differential Equations".
880 .. [4] `Cauchy-Riemann equations
881 <https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
882 Wikipedia.
884 Examples
885 --------
886 In the first example, we solve Bratu's problem::
888 y'' + k * exp(y) = 0
889 y(0) = y(1) = 0
891 for k = 1.
893 We rewrite the equation as a first-order system and implement its
894 right-hand side evaluation::
896 y1' = y2
897 y2' = -exp(y1)
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)