hexed::math Namespace Reference

Miscellaneous mathematical functions that aren't in std::math More...

Classes
class	Approx_equal
	functor to compare whether values are approximately equal More...
struct	Ball
	minimal representation of an n_dim-dimensional ball More...
struct	Root_options
	provides a convenient way to pass options to root-finding algorithms More...

Functions
template<typename number_t >
constexpr number_t	pow (number_t base, int exponent)
	Raises an arbitrary arithmetic type to an integer (not necessarily positive) power.
constexpr int	log (int base, int arg)
	Integer logarithm.
constexpr int	sign (bool condition)
	returns 1 if condition is true, otherwise -1
Eigen::VectorXi	direction (int n_dim, int i_dim, bool is_positive)
	the unit vector describing the direction from the center of an n_dim-dimensional Cartesian element to the face described by i_dim and sign
Eigen::VectorXi	direction (int n_dim, int i_face)
	the unit vector describing the direction from the center of an n_dim-dimensional Cartesian element to the i_faceth face
double	broyden (std::function< double(double)> error, double init_guess, Root_options opts, double init_diff=1e-3)
	Finds a root of a scalar function with Broyden's method.
double	bisection (std::function< double(double)> error, std::array< double, 2 > bounds, Root_options opts)
	Finds a root of a scalar function with the bisection method.
Mat	newton (std::function< Mat< dyn, dyn >(Mat<>)> error_jacobian, Mat<> guess, Root_options options)
	Finds a root of a vector function with Newton's method.
Eigen::VectorXd	hypercube_matvec (const Eigen::MatrixXd &, const Eigen::VectorXd &)
	Multiply every dimension of a (flattened) N-dimensional array by a matrix.
Eigen::VectorXd	dimension_matvec (const Eigen::MatrixXd &, const Eigen::VectorXd &, int i_dim)
	Multiply a single dimension of a hypercubic ND array by a matrix.
Eigen::VectorXd	pow_outer (const Eigen::VectorXd &, int n_dim)
	Raises a vector to a power via ND outer products.
template<int n_dim>
Eigen::Matrix< double, n_dim, n_dim >	orthonormal (Eigen::Matrix< double, n_dim, n_dim > basis, int i_dim)
	Orthonormalize a vector basis (with dimension \(\le 3\)).
Eigen::MatrixXd	orthonormal (Eigen::MatrixXd basis, int i_dim)
int	stretched_ind (int n_dim, int ind, std::array< bool, 2 > stretch)
	for indexing faces/vertices in Refined_face s with possible stretching
template<int n_dim>
double	interp (Mat< pow(2, n_dim)> values, Mat< n_dim > coords)
	\(n\)-linear interpolation of values.
template<typename vec_t >
vec_t	proj_to_segment (std::array< vec_t, 2 > endpoints, vec_t target)
template<typename T >
Mat	to_mat (T begin, T end)
	Constructs an Eigen::VectorXd from iterators begin() and end() to arithmetic types.
template<typename T >
Mat	to_mat (const T &range)
	Constructs an Eigen::VectorXd from any object supporting begin() and end() members.
template<int n_dim, int n_point>
Ball< n_dim >	bounding_ball (Mat< n_dim, n_point > points)
	computes a bounding ball of the convex hull of a set of points
template<int n_dim>
bool	intersects (Ball< n_dim > b, Mat< n_dim > endpoint0, Mat< n_dim > endpoint1)
	returns true if the ball b intersects the line through endpoint0 and endpoint1
double	chebyshev_step (int n_steps, int i_step, double safety=.9)
	the scaling factor for Chebyshev polynomial multistge time stepping
std::vector< double >	correct_values (std::vector< double > estimates, std::vector< double > exacts, double tol=huge)

Detailed Description

Miscellaneous mathematical functions that aren't in std::math

Function Documentation

bisection()

double hexed::math::bisection	(	std::function< double(double)>	error,
		std::array< double, 2 >	bounds,
		Root_options	opts )

Finds a root of a scalar function with the bisection method.

This is slower than broyden() but very robust.

Parameters

error	function to find the root of
bounds	Lower and upper bounds for a root.
opts	Use this to set the termination condition

bounding_ball()

template<int n_dim, int n_point>

Ball< n_dim > hexed::math::bounding_ball

(

Mat< n_dim, n_point >

points

)

computes a bounding ball of the convex hull of a set of points

Returns the smallest bounding ball centered at the center of mass of the points. This bounding ball is non-optimal but quick to compute. Note that a simplex is the convex hull of its vertices.

broyden()

double hexed::math::broyden	(	std::function< double(double)>	error,
		double	init_guess,
		Root_options	opts,
		double	init_diff = 1e-3 )

Finds a root of a scalar function with Broyden's method.

Parameters

error	error function to find the root of
init_guess	Initial guess for the root.
opts	Use this to set the termination condition
init_diff	How far away the second point used to initialize the derivative estimate should be.

chebyshev_step()

double hexed::math::chebyshev_step	(	int	n_steps,
		int	i_step,
		double	safety = .9 )

the scaling factor for Chebyshev polynomial multistge time stepping

Returns the ratio of the i_stepth step of n_steps to the nominal time step. A value of safety = 1 gives you the maximum stable time step and safety = 0 gives you a time step of 0 (although it isn't a linear function).

correct_values()

std::vector< double > hexed::math::correct_values	(	std::vector< double >	estimates,
		std::vector< double >	exacts,
		double	tol = huge )

Suppose that you compute some values, estimates, with a method that is very robust but has low accuracy. Suppose you also recompute these values as exacts with a method that is very accurate but not robust—It may give you extraneous values and/or fail to produce some of the correct values. The purpose of this function is to filter out the implausible values from exacts and give you only the ones that appear to be correct. The return value will have the same number of values as estimates, but some of them will be replaced with values from exacts in a way that obtains the best possible agreement. Replacements will be considered only if the difference between the "exact" value and the "estimate" is less than tol. The order of exacts does not matter. Values will be returned in the same order as estimates.

dimension_matvec()

Eigen::VectorXd hexed::math::dimension_matvec	(	const Eigen::MatrixXd &	mat,
		const Eigen::VectorXd &	vec,
		int	i_dim )

Multiply a single dimension of a hypercubic ND array by a matrix.

C.f. hypercube_matvec. If matrix is square, the shape of the output array will match the input. If matrix is a row vector, the output will still be hypercubic, but with one less dimension than the input. Otherwise, the resulting array will no longer be hypercubic.

hypercube_matvec()

Eigen::VectorXd hexed::math::hypercube_matvec	(	const Eigen::MatrixXd &	mat,
		const Eigen::VectorXd &	vec )

Multiply every dimension of a (flattened) N-dimensional array by a matrix.

Size of array along each dimension must be equal (i.e. the array is hypercube-shaped, or in my terminology, "hypercubic"). Matrix does not have to be square. Dimensionality of the array is inferred automatically by comparing the number of matrix columns to the array size. Array size must be a power of the number of matrix columns.

interp()

template<int n_dim>

double hexed::math::interp	(	Mat< pow(2, n_dim)>	values,
		Mat< n_dim >	coords )

\(n\)-linear interpolation of values.

ND generalization of bilinear interpolation.

Parameters

values	Values to interpolate. Assumed to be at corners of the unit hypercube.
coords	Coordinates to interpolate to.

log()

int hexed::math::log ( int base, int arg )

constexpr

Integer logarithm.

If base < 2, returns -1 to indicate failure. Otherwise, if arg < 1, returns 0. In the usual case where neither of the above are true, returns \(\lceil\log_{\mathtt{base}}(\mathtt{arg})\rceil\).

newton()

Mat hexed::math::newton	(	std::function< Mat< dyn, dyn >(Mat<>)>	error_jacobian,
		Mat<>	guess,
		Root_options	options )

Finds a root of a vector function with Newton's method.

Parameters

error_jacobian	Function to find the root of. If the input ( \( x \)) has \( n \) entries, this should be an \( n \times n + 1 \) matrix. The first column should be the error vector, and the remaing columns should be the Jacobian matrix, such that entry i, j of the Jacobian (i, j + 1 of the return matrix) should be the derivative of the ith component of the error w.r.t. the jth component of the input.
guess	initial guess for the root
options	use this to set the termination condition

orthonormal()

template<int n_dim>

Eigen::Matrix< double, n_dim, n_dim > hexed::math::orthonormal	(	Eigen::Matrix< double, n_dim, n_dim >	basis,
		int	i_dim )

Orthonormalize a vector basis (with dimension \(\le 3\)).

Assumes basis is invertible. Returns a matrix with the following properties:

• Unitary.
• Span of columns excluding the i_dimth is the same as for basis.
• Minimizes RMS difference between columns excluding i_dimth of return matrix and basis (sensitive to order).
• Inner product of i_dimth columns of return matrix and basis is positive.

pow()

template<typename number_t >

number_t hexed::math::pow ( number_t base, int exponent )

constexpr

Raises an arbitrary arithmetic type to an integer (not necessarily positive) power.

Can return constexpr, which std::pow is not allowed to do according to the standard (although the GCC implementation can anyway).

pow_outer()

Eigen::VectorXd hexed::math::pow_outer	(	const Eigen::VectorXd &	vec,
		int	n_dim )

Raises a vector to a power via ND outer products.

That is, takes an outer product with the vector {1} n_dim times along different dimensions to produce an n_dim-dimensional hypercubic array.

proj_to_segment()

template<typename vec_t >

vec_t hexed::math::proj_to_segment	(	std::array< vec_t, 2 >	endpoints,
		vec_t	target )

Finds the nearest point to target on the line segment defined by endpoints. Works for 2D or 3D.