mappy_ds.continuation.arclength_continuation

mappy_ds.continuation.arclength_continuation(fun, u0, *, step=0.01, step_min=1e-06, step_max=0.1, max_steps=1000, direction=1, bounds=None, newton_tol=1e-10, newton_max_iter=10, callback=None)

Follow a solution curve of F(u) = 0 by pseudo-arclength continuation.

F : R^{m+1} → R^m has a one-dimensional zero set. Starting from u0 (which is first corrected onto the curve with its last component held fixed), each step predicts u = u_k + h t_k along the unit tangent t_k and corrects it with Newton’s method on the bordered system

\[F(u) = 0, \qquad t_k \cdot (u - u_k) - h = 0,\]

whose Jacobian [J; t_kᵀ] is regular even where J loses rank in the parameter direction, so folds (turning points) are passed without special treatment (Keller, 1977). The step h grows by 1.5 when Newton converges in at most three iterations, is halved when six or more are needed or when it fails, and the orientation is kept by requiring t_k · t_{k+1} > 0.

Parameters:

funcallable

fun(u) -> (F, J) returning the residual (m,) and its Jacobian (m, m + 1) at u. It may raise MappyError or numpy.linalg.LinAlgError to signal that u is not evaluable; the step is then retried with a smaller h.

u0array_like

Starting point, shape (m + 1,); its last component is the parameter-like coordinate used by direction and bounds.

stepfloat, optional

Initial arclength step (Euclidean distance in the u space), by default 0.01.

step_min, step_maxfloat, optional

Bounds of the adaptive step, by default 1e-6 and 0.1. The continuation stops when a step below step_min would be needed.

max_stepsint, optional

Maximum number of accepted points after the first, by default 1000.

directionint, optional

+1 to start with the last component of u increasing, -1 decreasing, by default 1.

boundstuple of float, optional

Stop when the last component of u would leave [lo, hi]; the step is shortened so that the last accepted point lies close to the boundary.

newton_tolfloat, optional

Convergence tolerance of the corrector (on the Newton update relative to 1 + |u| and on |F|), by default 1e-10.

newton_max_iterint, optional

Maximum Newton iterations per step, by default 10.

callbackcallable, optional

Called with each accepted ArclengthPoint.

Returns:

list of ArclengthPoint

The accepted points in order, starting with the corrected u0. The list ends when bounds are left, max_steps is reached, the step cannot be reduced further, or the curve closes on itself.

Notes

The arclength is measured in the Euclidean metric of u; if the state and the parameter have very different scales, rescale the parameter in fun (or choose step accordingly) so that progress along the parameter is not starved.

Examples

The unit circle u0² + u1² = 1 traced from (1, 0) with u1 increasing first; it is closed, so the continuation stops after one turn.

>>> import numpy as np
>>> from mappy_ds.continuation import arclength_continuation
>>> def circle(u):
...     return np.array([u @ u - 1.0]), 2 * u[np.newaxis, :]
>>> pts = arclength_continuation(circle, [1.0, 0.0], step=0.3, step_max=0.3)
>>> all(abs(p.u @ p.u - 1) < 1e-9 for p in pts)
True
>>> pts[1].u[1] > 0 and 15 < len(pts) < 40
True