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.