NumCalc Namespace Reference

Detailed Description

Numerical integration and differentiation routines.

Classes
struct	QintCoefs
	Quadrature integration coefficients for a given number of points N. cq holds the weights; the step contribution is multiplied by dq_inv. Specialisations are provided for N = 1, 3, 5, 7, 9, 11, 13. More...

Enumerations
enum	Direction { zero_to_r , r_to_inf }
	Integration direction for additivePIntegral. More...
enum	t_grid { linear , logarithmic }
	Grid type for num_integrate. More...

Functions
template<typename C , typename... Args>
double	integrate (const double dt, std::size_t beg, std::size_t end, const C &f1, const Args &...rest)
	Quadrature integration of one or more vectors over a regular grid in t.
template<typename T >
std::vector< T >	derivative (const std::vector< T > &f, const std::vector< T > &drdt, const T dt, const int order=1)
	Computes the derivative df/dr on a non-uniform grid.
template<Direction direction, typename Real >
void	additivePIntegral (std::vector< Real > &answer, const std::vector< Real > &f, const std::vector< Real > &g, const std::vector< Real > &h, const Grid &gr, std::size_t pinf=0)
	Additively accumulates a partial integral into answer.
template<Direction direction, typename Real >
void	additivePIntegral (std::vector< Real > &answer, const std::vector< Real > &f, const std::vector< Real > &g, const Grid &gr, std::size_t pinf=0)
	Additively accumulates a partial integral into answer (two-function overload).
template<Direction direction, typename Real >
std::vector< Real >	partialIntegral (const std::vector< Real > &f, const std::vector< Real > &g, const std::vector< Real > &h, const Grid &gr, std::size_t pinf=0)
	Returns the partial integral ‘f(r) * Int[g(r’)*h(r'), {r', 0, r}]`.
std::function< double(long unsigned)>	linx (double a, double dt)
	Returns a function giving linearly-spaced x values: x(i) = a + i*dt.
std::function< double(long unsigned)>	one ()
	Returns a function that always returns 1.0 (uniform Jacobian)
std::function< double(long unsigned)>	logx (double a, double dt)
	Returns a function giving logarithmically-spaced x values: x(i) = a*exp(i*dt)
double	num_integrate (const std::function< double(double)> &f, double a, double b, long unsigned n_pts, t_grid type=linear)
	Numerically integrates a function f(x) over [a, b].

Variables
constexpr std::size_t	Nquad = 13
	Quadrature integration order; must be odd and in [1, 13].
constexpr QintCoefs< Nquad >	quintcoef
constexpr auto	cq = quintcoef.cq
constexpr auto	dq_inv = quintcoef.dq_inv

---

Class Documentation

NumCalc::QintCoefs

struct NumCalc::QintCoefs

Enumeration Type Documentation

Direction

enum NumCalc::Direction

Integration direction for additivePIntegral.

Enumerator
zero_to_r	Integrate from 0 outward to r.
r_to_inf	Integrate from infinity inward to r.

t_grid

enum NumCalc::t_grid

Grid type for num_integrate.

Function Documentation

integrate()

template<typename C , typename... Args>

double NumCalc::integrate ( const double  dt, std::size_t  beg, std::size_t  end, const C &  f1, const Args &...  rest  )

inline

Quadrature integration of one or more vectors over a regular grid in t.

Integrates the element-wise product f1[i] * rest[i] * ... from index beg to end-1 (i.e., end is exclusive), multiplied by the step size dt.

The grid must be uniform in the parametric variable t, but not necessarily in the physical variable x. For a non-uniform x grid, fold the Jacobian (dx/dt) into the integrand: pass f(x) * (dx/dt) rather than f(x) alone.

End-point corrections (using Nquad-point quadrature weights) are applied at the start and end of the range. The additivity property is preserved: int(a->b) + int(b->c) == int(a->c).

No safety checks are performed. Requirements:

• end - beg > 2 * Nquad
• end >= Nquad
• beg + 2 * Nquad <= f1.size()

Parameters

dt	Step size in the parametric variable t.
beg	First grid index (inclusive).
end	Last grid index (exclusive); if 0, defaults to f1.size().
f1	First vector to integrate (include Jacobian if grid is non-uniform in x).
rest	Additional vectors; all multiplied element-wise with f1.

Returns

Integral of the element-wise product times dt.

derivative()

template<typename T >

std::vector< T > NumCalc::derivative ( const std::vector< T > &  f, const std::vector< T > &  drdt, const T  dt, const int  order = 1  )

inline

Computes the derivative df/dr on a non-uniform grid.

Given f sampled on a non-uniform grid r(t) with uniform step dt, computes df/dr using finite difference coefficients:

df/dr = (df/di) / (dt * dr/dt)

where i is the grid index. Coefficients from: http://en.wikipedia.org/wiki/Finite_difference_coefficient

Uses a 7-point stencil in the interior; lower-order one-sided differences near the endpoints.

For order > 1, applies the derivative recursively.

Parameters

f	Function values on the grid.
drdt	Jacobian dr/dt at each grid point.
dt	Uniform step size in the parametric variable t.
order	Derivative order (default 1); higher orders applied recursively.

Returns

df/dr at each grid point.

additivePIntegral() [1/2]

template<Direction direction, typename Real >

void NumCalc::additivePIntegral ( std::vector< Real > &  answer, const std::vector< Real > &  f, const std::vector< Real > &  g, const std::vector< Real > &  h, const Grid &  gr, std::size_t  pinf = 0  )

inline

Additively accumulates a partial integral into answer.

Computes and adds to answer:

‘answer(r) += f(r) * Int[g(r’) * h(r'), {r', 0, r}]`

(or from infinity, depending on direction). Uses the trapezoidal rule.

Note

This is additive (+=). answer must already exist and be initialised (typically to zero).

Parameters

answer	Accumulation vector (modified in place).
f	Prefactor at each grid point.
g	First integrand factor.
h	Second integrand factor.
gr	Radial grid (provides dr/du and du).
pinf	Upper index limit; 0 means use full grid.

additivePIntegral() [2/2]

template<Direction direction, typename Real >

void NumCalc::additivePIntegral ( std::vector< Real > &  answer, const std::vector< Real > &  f, const std::vector< Real > &  g, const Grid &  gr, std::size_t  pinf = 0  )

inline

Additively accumulates a partial integral into answer (two-function overload).

As above, but integrand has only one factor g (no h):

‘answer(r) += f(r) * Int[g(r’), {r', 0, r}]`

Note

This is additive (+=). answer must already exist and be initialised (typically to zero).

partialIntegral()

template<Direction direction, typename Real >

std::vector< Real > NumCalc::partialIntegral	(	const std::vector< Real > &	f,
		const std::vector< Real > &	g,
		const std::vector< Real > &	h,
		const Grid &	gr,
		std::size_t	pinf = 0
	)

Returns the partial integral ‘f(r) * Int[g(r’)*h(r'), {r', 0, r}]`.

Allocates and returns the result vector; wrapper around additivePIntegral().

linx()

std::function< double(long unsigned)> NumCalc::linx ( double  a, double  dt  )

inline

Returns a function giving linearly-spaced x values: x(i) = a + i*dt.

one()

std::function< double(long unsigned)> NumCalc::one ( )

inline

Returns a function that always returns 1.0 (uniform Jacobian)

logx()

std::function< double(long unsigned)> NumCalc::logx ( double  a, double  dt  )

inline

Returns a function giving logarithmically-spaced x values: x(i) = a*exp(i*dt)

num_integrate()

double NumCalc::num_integrate ( const std::function< double(double)> &  f, double  a, double  b, long unsigned  n_pts, t_grid  type = linear  )

inline

Numerically integrates a function f(x) over [a, b].

Uses the same quadrature scheme as integrate(), on either a linear or logarithmic grid of n_pts points.

Returns 0 if a >= b or n_pts <= 1.

Parameters

f	Function to integrate.
a	Lower bound.
b	Upper bound.
n_pts	Number of quadrature points.
type	Grid spacing: linear (default) or logarithmic.

Returns

Approximate integral of f over [a, b].

Variable Documentation

Nquad

constexpr std::size_t NumCalc::Nquad = 13

constexpr

Quadrature integration order; must be odd and in [1, 13].