To provide a uniformed framework to root finding functions, we use several structures for storing different kind of functions. The pointer params is used to store the extra parameters. These new types come with dedicated macros starting in PNL_EVAL to evaluate the function and their Jacobian.
/* * f: R --> R * The function pointer returns f(x) * typedef struct { double (*F) (double x, void *params); void *params; } PnlFunc ; #define PNL_EVAL_FUNC(Fstruct, x) (*((Fstruct)->F))(x, (Fstruct)->params)
/* * f: R^2 --> R * The function pointer returns f(x) * typedef struct { double (*F) (double x, double y, void *params); void *params; } PnlFunc2D ; #define PNL_EVAL_FUNC2D(Fstruct, x, y) (*((Fstruct)->F))(x, y, (Fstruct)->params)
/* * f: R --> R * The function pointer computes f(x) and Df(x) and stores them in fx * and dfx respectively * typedef struct { void (*F) (double x, double *fx, double *dfx, void *params); void *params; } PnlFuncDFunc ; #define PNL_EVAL_FUNC_FDF(Fstruct, x, fx, dfx) (*((Fstruct)->F))(x, fx, dfx, (Fstruct)->params)
/* * f: R^n --> R * The function pointer returns f(x) * typedef struct { double (*F) (const PnlVect *x, void *params); void *params; } PnlRnFuncR ; #define PNL_EVAL_RNFUNCR(Fstruct, x) (*((Fstruct)->F))(x, (Fstruct)->params)
/* * f: R^n --> R^m * The function pointer computes the vector f(x) and stores it in * fx (vector of size m) * typedef struct { void (*F) (const PnlVect *x, PnlVect *fx, void *params); void *params; } PnlRnFuncRm ; #define PNL_EVAL_RNFUNCRM(Fstruct, x, fx) (*((Fstruct)->F))(x, fx, (Fstruct)->params) /* * Synonymous of PnlRnFuncRm for f:R^n --> R^n * typedef PnlRnFuncRm PnlRnFuncRn; #define PNL_EVAL_RNFUNCRN PNL_EVAL_RNFUNCRM
/* * f: R^n --> R^m * The function pointer computes the vector f(x) and stores it in fx * (vector of size m) * The Dfunction pointer computes the matrix Df(x) and stores it in dfx * (matrix of size m x n) * typedef struct { void (*F) (const PnlVect *x, PnlVect *fx, void *params); void (*DF) (const PnlVect *x, PnlMat *dfx, void *params); void (*FDF) (const PnlVect *x, PnlVect *fx, PnlMat *dfx, void *params); void *params; } PnlRnFuncRmDFunc ; #define PNL_EVAL_RNFUNCRM_DF(Fstruct, x, dfx) \ (*((Fstruct)->Dfunction))(x, dfx, (Fstruct)->params) #define PNL_EVAL_RNFUNCRM_FDF(Fstruct, x, fx, dfx) \ (*((Fstruct)->F))(x, fx, dfx, (Fstruct)->params) #define PNL_EVAL_RNFUNCRM_F_DF(Fstruct, x, fx, dfx) \ if ( (Fstruct)->FDF != NULL ) \ { \ PNL_EVAL_RNFUNCRN_FDF (Fstruct, x, fx, dfx); \ } \ else \ { \ PNL_EVAL_RNFUNCRN (Fstruct, x, fx); \ PNL_EVAL_RNFUNCRN_DF (Fstruct, x, dfx); \ } /* * Synonymous of PnlRnFuncRmDFunc for f:R^n --> R^m * typedef PnlRnFuncRmDFunc PnlRnFuncRnDFunc; #define PNL_EVAL_RNFUNCRN_DF PNL_EVAL_RNFUNCRM_DF #define PNL_EVAL_RNFUNCRN_FDF PNL_EVAL_RNFUNCRM_FDF #define PNL_EVAL_RNFUNCRN_F_DF PNL_EVAL_RNFUNCRM_F_DF
To use the following functions, you should include pnl/pnl_root.h.
Real valued functions of a real argument
double pnl_root_brent (PnlFunc *F, double x1, double x2, double *tol)
Description Find the root of F between x1 and x2 with an accuracy of order tol. On
exit tol is an upper bound of the error.
int pnl_root_newton_bisection (PnlFuncDFunc *Func, double x_min, double
x_max, double tol, int N_Max, double *res)
Description Find the root of F between x1 and x2 with an accuracy of order tol
and a maximum of N_max iterations. On exit, the root is stored in res. Note that the
function Func must also compute the first derivative of the function. This function relies
on combining Newton’s approach with a bisection technique.
int pnl_root_newton (PnlFuncDFunc *Func, double x0, double x_eps, double fx_eps, int
max_iter, double *res)
Description Find the root of f starting from x0 using Newton’s method with descent
direction given by the inverse of the derivative, ie. dk = f(xk)∕f′(xk). Armijo’s line search is
used to make sure |f| decreases along the iterations. αk = max{γj ; j ≥ 0} such
that
| |f(xk + αkdk)| | ≤|f(xk)|(1 -ωαk). |
In this implementation, ω = 10-4 and γ = 1∕2. The algorithm stops when one of the three following conditions is met:
the maximum number of iterations max_iter is reached;
the last improvement over x is smaller that x * x_eps;
at the current position |f(x)| < fx_eps
On exit, the root is stored in res.
int pnl_root_bisection (PnlFunc *Func, double xmin, double xmax, double epsrel, double
espabs, int N_max, double *res)
Description Find the root of F between x1 and x2 with the accuracy |x2 - x1| < epsrel * x1
+ epsabs or with the maximum number of iterations N_max. On exit, res = (x2 + x1) /
2.
Vector valued functions with several arguments
int pnl_multiroot_newton (PnlRnFuncRnDFunc *func, const PnlVect *x0, double
x_eps, double fx_eps, int max_iter, int verbose, PnlVect *res)
Description Find the root of func starting from x0 using Newton’s method with descent
direction given by the inverse of the Jacobian matrix, ie. dk = (∇f(xk))-1f(xk). Armijo’s line
search is used to make sure |f| decreases along the iterations. αk = max{γj ; j ≥ 0} such
that
| |f(xk + αkdk)| | ≤|f(xk)|(1 -ωαk). |
In this implementation, ω = 10-4 and γ = 1∕2. The algorithm stops when one of the three following conditions is met:
the maximum number of iterations max_iter is reached;
the norm of the last improvement over x is smaller that |x| * x_eps;
at the current position |f(x)| < fx_eps
On exit, the root is stored in res. Note that the function F must also compute the first derivative of the function. When defining Func, you must either define Func->F and Func->DF or Func->FDF.
We provide two wrappers for calling minpack routines.
int pnl_root_fsolve (PnlRnFuncRnDFunc *f, PnlVect *x, PnlVect *fx, double xtol, int
maxfev, int *nfev, PnlVect *scale, int error_msg)
Description Compute the root of a function f : ℝn⟼ℝn. Note that the number of
components of f must be equal to the number of variates of f. This function returns OK or
FAIL if something went wrong.
Parameters
f is a pointer to a PnlRnFuncRnDFunc used to store the function whose root is to be found. f can also store the Jacobian of the function, if not it is computed using finite differences (see the file examples/minpack_test.c for a usage example). f->FDF can be NULL because it is not used in this function.
x contains on input the starting point of the search and an approximation of the root of f on output,
xtol is the precision required on x, if set to 0 a default value is used.
maxfev is the maximum number of evaluations of the function f before the algorithm returns, if set to 0, a coherent number is determined internally.
nfev contains on output the number of evaluations of f during the algorithm,
scale is a vector used to rescale x in a way that each coordinate of the solution is approximately of order 1 after rescaling. If on input scale=NULL, a scaling vector is computed internally by the algorithm.
error_msg is a boolean (TRUE or FALSE) to specify if an error message should be printed when the algorithm stops before having converged.
On output, fx contains f(x).
int pnl_root_fsolve_lsq (PnlRnFuncRmDFunc *f, PnlVect *x, int m, PnlVect *fx,
double xtol, double ftol, double gtol, int maxfev, int *nfev, PnlVect *scale, int
error_msg)
Description Compute the root of x ∈ℝn⟼∑
i=1mfi(x)2, note that there is no reason why m
should be equal to n.
Parameters
f is a pointer to a PnlRnFuncRmDFunc used to store the function whose root is to be found. f can also store the Jacobian of the function, if not it is computed using finite differences (see the file examples/minpack_test.c for a usage example). f->FDF can be NULL because it is not used in this function.
x contains on input the starting point of the search and an approximation of the root of f on output,
m is the number of components of f,
xtol is the precision required on x, if set to 0 a default value is used.
ftol is the precision required on f, if set to 0 a default value is used.
gtol is the precision required on the Jacobian of f, if set to 0 a default value is used.
maxfev is the maximum number of evaluations of the function f before the algorithm returns, if set to 0, a coherent number is determined internally.
nfev contains on output the number of evaluations of f during the algorithm,
scale is a vector used to rescale x in a way that each coordinate of the solution is approximately of order 1 after rescaling. If on input scale=NULL, a scaling vector is computed internally by the algorithm.
error_msg is a boolean (TRUE or FALSE) to specify if an error message should be printed when the algorithm stops before having converged.
On output, fx contains f(x).