Math

Nonlinear solver

template<class F, class T = double>
class Solve

A class that solves a provided function for a zero.

The first template argument, F, is the type of the function. Typically you would simple function object that defines the operator() method and use that for F.

The second template argument (optional: default = double) is the type of the arguments to the function.

The solving process is in two parts:

1. First the solution needs to be bracketed.

   This can be done in several ways:

   a) If you know the appropriate range, you can set it in the constructor:

   Solve<MyFunc> solver(func, x_min, x_max);
   

   b) If you don’t know the range, you can set an initial guess range and let Solve expand the range until it finds a range that brackets the solution:

   Solve<MyFunc> solver(func, x_min, x_max);
   solver.bracket()
   

   c) If you know one end of the range, but not the other you can expand only one side of the range from the initial guess:

   Solve<MyFunc> solver(func, 0, r_max);
   solver.bracketUpper()
   

   (There is also a corresponding bracketLower() method.)

   d) Sometimes there is a fundamental limit past which you want to make sure that the range doesn’t go. e.g. For values that must be positive, you might want to make sure the range doesn’t extend past 0.

   For this case, there are the following two functions

   solver.bracketUpperWithLimit(upper_limit)
   solver.bracketLowerWithLimit(lower_limit)
   

   Finally, note that there is nothing in the code that enforces x_min < x_max. If the two limits bracket the code in reverse order, that’s ok.

2. The next step is to solve for the root within this range. There are currently two algorithms for doing that: Bisect and Brent. You can tell Solve which one to use with:

   solver.setMethod(Bisect)
   solver.setMethod(Brent)
   

   (In fact, the former is unnecessary, since Bisect is the default.)

   Then the method root() will solve for the root in that range.

Typical usage:

struct MyFunc {
    double operator()(double x) { return [...] }
};

MyFunc func;
Solve<MyFunc> solver(func, xmin, xmax);
solver.bracket();         // If necessary
solver.setMethod(Brent);  // If desired
double result = solver.root()

Public Functions

inline Solve(const F &func_, T lb_, T ub_)

Constructor taking the function to solve, and the range (if known).

inline void setMaxSteps(int m)

Set the maximum number of steps for the solving algorithm to use.

inline void setMethod(Method m_)

Set which method to use: Bisect or Brent.

inline void setXTolerance(T tol)

Set the tolerance to define when the root is close enough to 0. (abs(x) < tol)

inline T getXTolerance() const

Get the current tolerance.

inline void setBounds(T lb, T ub)

Set the bounds for the search to new values.

inline T getLowerBound() const

Get the current search bounds.

inline T getUpperBound() const

inline void evaluateBounds() const

Make sure the current bounds have corresponding flower, fupper.

inline void bracket()

Hunt for bracket, geometrically expanding range.

This version assumes that the root is to the side of the end point that is closer to 0. This will be true if the function is monotonic.

inline bool bracket1(T &a, T &b, T &fa, T &fb)

inline void bracketUpper()

Hunt for bracket, geometrically expanding range.

This one only expands to the right for when you know that the lower bound is definitely to the left of the root, but the upper bound might not bracket it.

inline void bracketLower()

Hunt for bracket, geometrically expanding range.

The opposite of bracketUpper &#8212; only expand to the left.

inline bool bracket1WithLimit(T &a, T &b, T &fa, T &fb, T &c)

inline void bracketUpperWithLimit(T upper_limit)

Hunt for upper bracket, with an upper limit to how far it can go.

inline void bracketLowerWithLimit(T lower_limit)

Hunt for lower bracket, with a lower limit to how far it can go.

inline T root() const

Find the root according the the method currently set.

inline T bisect() const

A simple bisection root-finder.

inline T zbrent() const

A more sophisticated root-finder using Brent’s algorithm.

Bessel and Related Functions

double galsim::math::cyl_bessel_j(double nu, double x)

double galsim::math::cyl_bessel_y(double nu, double x)

double galsim::math::cyl_bessel_k(double nu, double x)

double galsim::math::cyl_bessel_i(double nu, double x)

double galsim::math::j0(double x)

double galsim::math::j1(double x)

double galsim::math::getBesselRoot0(int s)

double galsim::math::getBesselRoot(double nu, int s)

Other mathematical functions

void galsim::math::sincos(double theta, double &sint, double &cost)

double galsim::math::gamma_p(double a, double x)

double galsim::math::sinc(double x)

double galsim::math::Si(double x)

double galsim::math::Ci(double x)

Horner’s method for polynomial evaluation

void galsim::math::Horner(const double *x, const int nx, const double *coef, const int nc, double *result)

void galsim::math::Horner2D(const double *x, const double *y, const int nx, const double *coef, const int ncx, const int ncy, double *result, double *temp)

C++ Integration Functions

template<class T>
struct IntRegion

A type that encapsulates everything known about the integral in a region.

The constructor of IntRegion takes the minimum and maximum values for the region, along with an optional ostream for outputing diagnostic information.

After the integration is done, the IntRegion will also hold the estimate of the integral’s value over the region, along with an estimate of the error.

Public Functions

inline IntRegion(const T a, const T b, std::ostream *dbgout_ = 0, std::map<T, T> *fxmap_ = 0)

Constructor taking the bounds of the region and (optionally) an ostream for outputing diagnostic information.

Note that normally a < b. However, a > b is also valid. If a > b, then the integral will be from right to left, which is just the negative of the integral from b to a.

Parameters:

• a – The left end of the region

• b – The right end of the region

• dbgout_ – An optional ostream for diagnostic info

• fxmap_ – Known results

inline bool operator<(const IntRegion<T> &r2) const

op< sorts by the error estimate

inline bool operator>(const IntRegion<T> &r2) const

op> sorts by the error estimate

inline void subDivide(std::vector<IntRegion<T>> &children)

Subdivide a region according to the current split points or biset.

If there are no split points set yet, then this will bisect the region. However, you may also set split points using addSplit(x). This is worth doing if you know about any discontinuities, zeros or poles in the function you are integrating.

inline void bisect()

Set a split point at the bisection.

inline void findZeroCrossings()

Search for zero crossings based on values stored in fxmap.

inline void addSplit(const T x)

Add a split point to the current list to be used by the next subDivide call.

This is worth doing if you know about any discontinuities, zeros or poles in the function you are integrating.

inline size_t getNSplit() const

Get the number of split points currently set.

inline const T &left() const

Get the left end of the region.

inline const T &right() const

Get the right end of the region.

inline const T &getErr() const

Get the current error estimate.

inline const T &getArea() const

Get the current estimate of the integral over the region.

inline void setArea(const T &a, const T &e)

Set a new estimate of the area and error.

inline void useFXMap()

Setup an fxmap for this region.

Public Members

std::ostream *dbgout

std::map<T, T> *fxmap

template<class UF>
inline double galsim::integ::int1d(const UF &func, double min, double max, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 1-dimensional integral using simple min/max values for the region.

Parameters:

• func – The function to be integrated (may be a function object)

• min – The lower bound of the integration

• max – The upper bound of the integration

• relerr – The target relative error

• abserr – The target absolute error

template<class UF>
inline double galsim::integ::int1d(const UF &func, IntRegion<double> &reg, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 1-dimensional integral using an IntRegion.

Parameters:

• func – The function to be integrated(may be a function object)

• reg – The region of the integration

• relerr – The target relative error

• abserr – The target absolute error

template<class BF>
inline double galsim::integ::int2d(const BF &func, double xmin, double xmax, double ymin, double ymax, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 2-dimensional integral using simple min/max values for borh regions (i.e. the integral is over a square)

template<class BF, class YREG>
inline double galsim::integ::int2d(const BF &func, IntRegion<double> &reg, const YREG &yreg, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 2-dimensional integral.

Parameters:

• func – The function to be integrated (may be a function object)

• reg – The region of the inner integral

• yreg – yreg(x) is the region of the outer integral at x

• relerr – The target relative error

• abserr – The target absolute error

template<class BF>
inline double galsim::integ::int2d(const BF &func, IntRegion<double> &reg, IntRegion<double> &yreg, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 3-dimensional integral using constant IntRegions for both regions (i.e. the integral is over a square)

template<class TF>
inline double galsim::integ::int3d(const TF &func, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 3-dimensional integral using simple min/max values for all regions (i.e. the integral is over a cube)

template<class TF, class YREG, class ZREG>
inline double galsim::integ::int3d(const TF &func, IntRegion<double> &reg, const YREG &yreg, const ZREG &zreg, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 3-dimensional integral.

Parameters:

• func – The function to be integrated (may be a function object)

• reg – The region of the inner integral

• yreg – yreg(x) is the region of the middle integral at x

• zreg – zreg(x)(y) is the region of the outer integral at x,y

• relerr – The target relative error

• abserr – The target absolute error

template<class TF>
inline double galsim::integ::int3d(const TF &func, IntRegion<double> &reg, IntRegion<double> &yreg, IntRegion<double> &zreg, const double relerr = DEFRELERR, const double abserr = DEFABSERR)

Perform a 3-dimensional integral using constant IntRegions for all regions (i.e. the integral is over a cube)

double galsim::math::hankel_trunc(const std::function<double(double)> f, double k, double nu, double maxr, double relerr = 1.e-6, double abserr = 1.e-12, int nzeros = 10)

double galsim::math::hankel_inf(const std::function<double(double)> f, double k, double nu, double relerr = 1.e-6, double abserr = 1.e-12, int nzeros = 10)

Misc Utilities

template<typename T>
bool galsim::math::isNan(T x)

int galsim::SetOMPThreads(int num_threads)

int galsim::GetOMPThreads()