nutopy.path.solve

solve(f, x0, pars0, parsf, *, df=None, options=None, callback=None, args=(), kwargs={})[source]

Solver of one-parameter family of nonlinear equations:

\[f(x, \lambda) = 0, \quad \lambda \in [\lambda_0, \lambda_f]\]

with \(f : \mathbf{R}^n \times \mathbf{R} \rightarrow \mathbf{R}^n\).

The solver needs an initial pair \((x_0, \lambda_0)\) such that \(f(x_0, \lambda_0)=0\), that is the pair is a zero of \(f\). The goal is to compute a path of zeros \((x_n, \lambda_n)\) such that the final pair denoted \((x_N, \lambda_N)\) satisfies \(f(x_N, \lambda_N)=0\) and \(\lambda_N=\lambda_f\). The final point \(x_N\) is then denoted \(x_f\).

The argument \(\lambda\) is refered as a parameter. The solver admits a vector of parameters denoted pars instead of a scalar parameter. In that case, the solver solves

\[f(x, (1-\lambda)\, \texttt{pars}_0 + \lambda\, \texttt{pars}_f) = 0, \quad \lambda \in [0, 1],\]

and the goal is to compute a path of zeros whose extremities satisfy \(f(x_0, \texttt{pars}_0) = 0\) and \(f(x_f, \texttt{pars}_f) = 0\).

Notes

The differential path following algorithm of the nutopy.path.solve method is partially described in Ref 1. The method is a Predictor-Corrector method 2 with high-order Runge-Kutta schemes for the prediction and a simplified Newton for the correction. Only the prediction is descibred in 1.

Parameters
  • f (callable) – A vector function representing f(x, pars).

  • x0 (float vector or float) – Initial solution, that is f(x0, pars0) = 0.

  • pars0 (float vector or float) – Initial parameters.

  • parsf (float vector or float) – Final parameters, the goal is to compute xf such that F(xf, parsf) = 0.

  • df (callable, optional, key only) –

    The Jacobian f'(x, pars) which returns a tuple of float matrices. It must return a tuple (dfdx, dfdpars) where dfdx = df/dx(x, pars) and dfdpars = df/dpars(x,pars).

    dfdx has shape (n,n) where n is the size of x and dfdpars has shape (n,m) where m is the size of pars.

  • options (nutopy.path.Options, optional, key only) – A dictionnary of solver options.

  • callback (callable, optional, key only) –

    A function of the form callback(infos) which returns a flag (named status). The path.solve method stops if status < 0 at some iteration. The argument infos is a structure containing:

    • infos.x : x at current iteration;

    • infos.pars : pars at current iteration;

    • infos.f : f(x,pars) at current iteration.

    Advice. Return a status lower than -10 to avoid any confusion with status below.

  • args (tuple, key only) – optional arguments for f (and df if provided).

  • kwargs (dictionnary, key only) – keywords optional arguments for f (and df if provided).

Returns

sol – A dictionary/struct of outputs:

  • sol.xf or sol.get('xf') (float vector) : final solution xf;

  • sol.parsf (float vector) : final parameters parsf;

  • sol.xout (float matrix)solutions along the path of zeros,

    ie f(xout[k], parsout[k]) = 0;

  • sol.parsout (float matrix) : parameters at each computed step;

  • sol.fnorm (float vector) : norm of the function f at each computed step;

  • sol.status (int) : termination status of the solver

    • if status== 0 then MaxSf is too small;

    • if status==-1 then input is not consistent;

    • if status==-2 then larger MaxSteps is needed;

    • if status==-3 then step size became too small;

    • if status==-4 then problem is problably stiff;

    • if status==-5 then |f| is greater than NormeSfunMax;

    • if status==-6 then the homotopic parameter do not evolve relatively w.r.t. TolHomParEvolution;

    • if status==-7 then a turning point occurs.

  • sol.success (bool) : whether or not the solver exited successfully;

  • sol.message (str) : description of the cause of the termination;

  • sol.sout (float vector) : arc length at each computed step;

  • sol.vout (float matrix) : vout[:,k] = c'[:,k] where c = (x, pars);

  • sol.dets (float vector) : det(f'(c(s)),c'(s)) at each computed step;

  • sol.ps (float vector) :( c'(s_old) | c'(s) ) at each computed step.

Notes

  • Use numpy.ndarray to define float vectors.

  • You can use a Newton solver (see nutopy.nle.solve) to get x0.

  • You can define a unique function f coding f and its derivative:

>>> from nutopy.tools import *
>>> def df(x, dx, p, dp, a):
...     dy    = np.zeros(x.size)
...     dy[0] = a * dx[0] - dp[0]
...     dy[1] = np.exp(x[0] + x[1]) * ( dx[0] + dx[1] ) - dp[1]
...     return dy
>>> @tensorize(df, tvars=(1, 2))
... def f(x, p, a): # p is pars
...     y    = np.zeros(x.size)
...     y[0] = a * x[0] - p[0]
...     y[1] = np.exp(x[0] + x[1]) - p[1]
...     return y

Then, call path.solve with df=f:

>>> from nutopy import path
>>> x0 = np.array([0.0, 0.0])
>>> p0 = np.array([0.0, 1.0])
>>> pf = np.array([1.0, 2.0])
>>> a  = 2.0
>>> sol = path.solve(f, x0, p0, pf, args=a, df=f)
>>> # or equivalently: df is not required if f is tensorized
>>> sol = path.solve(f, x0, p0, pf, args=a)

Examples

>>> import numpy as np
>>> from nutopy import path

Example 1: scalar

>>> def f(x, l):
...     return 2.0*x-l
>>> x0  = 0.0
>>> l0  = 0.0
>>> lf  = 1.0
>>> sol = path.solve(f, x0, l0, lf)

Example 2: scalar, with numpy

>>> def f(x, l):
...     y    = np.zeros(x.size)
...     y[0] = 2.0*x[0] - l[0]
...     return y
>>> x0  = np.array([0.0])
>>> l0  = np.array([0.0])
>>> lf  = np.array([1.0])
>>> sol = path.solve(f, x0, l0, lf)

Example 3: with Jacobian, options and optional argument

>>> def f(x, l, a):
...     y    = np.zeros(x.size)
...     y[0] = a*x[0]-l[0]
...     y[1] = np.exp(x[0]+x[1])-1.0
...     return y
>>> def df(x, l, a):
...     n       = x.size
...     m       = l.size
...     dfdx    = np.zeros((n,n))
...     dfdl    = np.zeros((n,m))
...
...     dfdx[0,0] = a
...     dfdx[1,0] = np.exp(x[0]+x[1])
...     dfdx[1,1] = np.exp(x[0]+x[1])
...
...     dfdl[0,0] = -1.0
...     return (dfdx, dfdl)
>>> opt = path.Options(ODESolver='dopri5', MaxIterCorrection=0);
>>> a   = 2.0
>>> x0  = np.array([0.0, 0.0])
>>> l0  = 0.0
>>> lf  = 1.0
>>> sol = path.solve(f, x0, l0, lf, options=opt, args=a, df=df)

Example 4: homotopy on two parameters with Jacobian

>>> def f(x, l, a):
...     y = np.zeros(x.size)
...     y[0] = a * x[0] - l[0]
...     y[1] = np.exp(x[0] + x[1]) - l[1]
...     return y
>>> def df(x, l, a):
...     n = x.size
...     m = l.size
...     dfdx = np.zeros((n, n))
...     dfdl = np.zeros((n, m))
...
...     dfdx[0, 0] = a
...     dfdx[1, 0] = np.exp(x[0] + x[1])
...     dfdx[1, 1] = np.exp(x[0] + x[1])
...
...     dfdl[0, 0] = -1.0
...     dfdl[1, 1] = -1.0
...     return (dfdx, dfdl)
>>> x0  = np.array([0.0, 0.0])
>>> l0  = np.array([0.0, 1.0])
>>> lf  = np.array([1.0, 2.0])
>>> a   = 2.0
>>> sol = path.solve(f, x0, l0, lf, args=a, df=df)

Example 5: use callback to stop homotopy

>>> def f(x, l):
...     y = np.zeros(x.size)
...     y[0] = 2.0 * x[0] - l[0]
...     return y
>>> def mycallback(infos):
...     if(infos.x[0]>=1.0/4.0): # we stop if x >= 1/4
...         return -11
...     return 0
>>> x0 = np.array([0.0])
>>> l0 = np.array([0.0])
>>> lf = np.array([1.0])
>>> sol = path.solve(f, x0, l0, lf, callback=mycallback)
>>> print(sol.status)
-11

References

1(1,2)

J.-B. Caillau, O. Cots and J. Gergaud, Differential continuation for regular optimal control problems, Optimization Methods and Software, 27 (2011), no. 2, 177–196.

2

E. Allgower & K. Georg, Introduction to numerical continuation methods, vol. 45 of Classics in Applied Mathematics, Soc. for Industrial and Applied Math., Philadelphia, PA, USA, (2003), xxvi+388.

Raises