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 isf(x0, pars0) = 0
.pars0 (float vector or
float
) – Initial parameters.parsf (float vector or
float
) – Final parameters, the goal is to compute xf such thatF(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)
wheredfdx = df/dx(x, pars)
anddfdpars = df/dpars(x,pars)
.dfdx
has shape(n,n)
wheren
is the size ofx
anddfdpars
has shape(n,m)
wherem
is the size ofpars
.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 (namedstatus
). Thepath.solve
method stops ifstatus < 0
at some iteration. The argumentinfos
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
(anddf
if provided).kwargs (dictionnary, key only) – keywords optional arguments for
f
(anddf
if provided).
- Returns
sol – A dictionary/struct of outputs:
sol.xf
orsol.get('xf')
(float vector) : final solutionxf
;sol.parsf
(float vector) : final parametersparsf
;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 functionf
at each computed step;sol.status
(int
) : termination status of the solverif
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]
wherec = (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.
See also
Notes
Use
numpy.ndarray
to define float vectors.You can use a Newton solver (see
nutopy.nle.solve
) to getx0
.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
withdf=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
ArgumentTypeError – Variable x0 must be a one dimensional numpy array or a float.
ArgumentTypeError – Variable pars0 must be a one dimensional numpy array or a float.
ArgumentTypeError – Variable parsf must be a one dimensional numpy array or a float.
ArgumentTypeError – Variable options must be a Options object.
ArgumentDimensionError – Variable x0 must be a one dimensional numpy array or a float.
ArgumentDimensionError – Variable pars0 must be a one dimensional numpy array or a float.
ArgumentDimensionError – Variable parsf must be a one dimensional numpy array or a float.
ArgumentDimensionError – Variable parsf must have the same size than pars0.