Table of contents


Back to the main page Back to the documentation 1. Overview: GSL wrappers 2. ODESolver 3. Spline

Overview: GSL wrappers


We provide some simple wrappers on the GSL library for solving ODEs and making Splines. The classes take care of all the memory handling and provide a simple interface for performing the operations.

ODESolver


Solve coupled ODEs by providing the right hand side as a lambda (std::function):

using namespace FML::SOLVERS::ODESOLVER;

// The RHS of the ODE y'' + y = 0 written as y0' = y1 ; y1' = -y0
ODEFunction dydx = [&](double x, const double* y, double* dydx) {
  dydx[0] =  y[1];
  dydx[1] = -y[0];
  return GSL_SUCCESS;
};

// x-values we store the points at
DVector x_array = FML::MATH::linspace(xmin, xmax, npts);

// Initial value
DVector y_ini{1.0, 0.0};

// Solve the ODE over the points in x_array
ODESolver ode;
ode.solve(dydx, x_array, y_ini);

// You can also choose the method if you want
// ode.solve(dydx, x_array, y_ini, gsl_odeiv2_step_rk2);

// Get the result
auto y_array = ode.get_data_by_component(0);

Some methods require the Jacobian of the system to be provided, that is done in the same way. See the header file for more info.

Spline


Make a simple qubic (customizable) spline from arrays. Lookups are thread-safe as long as its not created inside a thread.

// Make a spline
DVector x_array = FML::MATH::linspace(xmin, xmax, npts);
DVector y_array = FML::MATH::linspace(xmin, xmax, npts);
Spline spline(x_array, y_array, "My Spline");

// Evaluate f(0)
double value = spline(0.0);

You can get warnings if you evaluate it out of bounds.

Make a 2D spline from arrays. Lookups are thread-safe as long as its not created inside a thread.

// Make a spline
DVector x_array   = FML::MATH::linspace(xmin, xmax, npts);
DVector y_array   = FML::MATH::linspace(xmin, xmax, npts);
DVector2D z_array(npts, npts);

// ... fill it ...

Spline2D spline(x_array, y_array, z_array, "My Spline");

// Evaluate f(0,0)
double value = spline(0.0, 0.0);

You can get warnings if you evaluate it out of bounds. You can also choose the interpolation method, see the header file and the GSL documentation for more info.