1 Using the C++ interface (oscode)¶

1.1 Overview¶

This documentation illustrates how one can use oscode via its C++ interface. Usage of oscode involves

defining an equation to solve,
solving the equation,
and extracting the solution and other statistics about the run.

The next sections will cover each of these. For a complete reference, see the C++ interface reference page, and for examples see the examples directory on GitHub.

1.2 Defining an equation¶

The equations oscode can be used to solve are of the form

\[\ddot{x}(t) + 2\gamma(t)\dot{x}(t) + \omega^2(t)x(t) = 0,\]

where \(x(t)\), \(\gamma(t)\), \(\omega(t)\) can be complex. We will call \(t\) the independent variable, \(x\) the dependent variable, \(\omega(t)\) the frequency term, and \(\gamma(t)\) the friction or first-derivative term.

Defining an equation is via

giving the frequency \(\omega(t)\),
giving the first-derivative term \(\gamma(t)\),

Defining the frequency and the first-derivative term can either be done by giving them as functions explicitly, or by giving them as sequences evaluated on a grid of \(t\).

1.2.1 \(\omega\) and \(\gamma\) as explicit functions¶

If \(\omega\) and \(\gamma\) are closed-form functions of time, then define them as

#include "solver.hpp" // de_system, Solution defined in here

std::complex<double> g(double t){
    return 0.0;
};

std::complex<double> w(double t){
    return std::pow(9999,0.5)/(1.0 + t*t);
};

Then feed them to the solver via the de_system class:

de_system sys(&w, &g);
Solution solution(sys, ...) // other arguments left out

1.2.2 \(\omega\) and \(\gamma\) as time series¶

Sometimes \(\omega\) and \(\gamma\) will be results of numerical integration, and they will have no closed-form functional form. In this case, they can be specified on a grid, and oscode will perform linear interpolation on the given grid to find their values at any timepoint. Because of this, some important things to note are:

oscode will assume the grid of timepoints \(\omega\) and \(\gamma\) are not evenly spaced. If the grids are evenly sampled, set even=true in the call for de_system(), this will speed linear interpolation up significantly.
The timepoints grid needs to be monotonically increasing.
The timepoints grid needs to include the range of integration (\(t_i\),:math:t_f).
The grids for the timepoints, frequencies, and first-derivative terms have to be the same size.
The speed/efficiency of the solver depends on how accurately it can carry out numerical integrals of the frequency and the first-derivative terms, therefore the grid fineness needs to be high enough. (Typically this means that linear interpolation gives a \(\omega(t)\) value that is accurate to 1 part in \(10^{6}\) or so.) If you want oscode to check whether the grids were sampled finely enough, set check_grid=true in the call for de_system().

To define the grids, use any array-like container which is contiguous in memory, e.g. an Eigen::Vector, std::array, std::vector:

#include "solver.hpp" // de_system, Solution defined in here

// Create a fine grid of timepoints and
// a grid of values for w, g
N = 10000;
std::vector<double> ts(N);
std::vector<std::complex<double>> ws(N), gs(N);

// Fill up the grids
for(int i=0; i<N; i++){
    ts[i] = i;
    ws[i] = std::sqrt(i);
    gs[i] = 0.0;
}

They can then be given to the solver again by feeding a pointer to their underlying data to the de_system class:

de_system sys(ts.data(), ws.data(), gs.data());
Solution solution(sys, ...) // other arguments left out

Often \(\omega\) and \(\gamma\) are much easier to perform linear interpolation on once taken natural log of. This is what the optional islogw and islogg arguments of the overloaded de_system::de_system() constructor are for:

#include "solver.hpp" // de_system, Solution defined in here

// Create a fine grid of timepoints and
// a grid of values for w, g
N = 10000;
std::vector<double> ts(N);
std::vector<std::complex<double> logws(N), gs(N); // Note the log!

// Fill up the grids
for(int i=0; i<N; i++){
    ts[i] = i;
    logws[i] = 0.5*i;
    gs[i] = 0.0; // Will not be logged
}

// We want to tell de_system that w has been taken natural log of, but g
// hasn't. Therefore islogw=true, islogg=false:
de_system sys(ts.data(), logws.data(), gs.data(), true, false);
Solution solution(sys, ... ) // other arguments left out

1.2.2.1 DIY interpolation¶

For some problems, linear interpolation of \(\omega\) and \(\gamma\) (or their natural logs) might simply not be enough.

For example, the user could carry out cubic spline interpolation and feed \(\omega\) and \(\gamma\) as functions to de_system.

Another example for wanting to do (linear) interpolation outside of oscode is when Solution.solve() is ran in a loop, and for each iteration a large grid of \(\omega\) and \(\gamma\) is required, depending on some parameter. Instead of generating them over and over again, one could define them as functions, making use of some underlying vectors that are independent of the parameter we iterate over:

// A, B, and C are large std::vectors, same for each run
// k is a parameter, different for each run
// the grid of timepoints w, g are defined on starts at tstart, and is
// evenly spaced with a spacing tinc.

// tstart, tinc, A, B, C defined here

std::complex<double> g(double t){
    int i;
    i=int((t-tstart)/tinc);
    std::complex<double> g0 = 0.5*(k*k*A[i] + 3.0 - B[i] + C[i]*k;
    std::complex<double> g1 = 0.5*(k*k*A[i+1] + 3.0 - B[i+1] + C[i+1]*k);
    return (g0+(g1-g0)*(t-tstart-tinc*i)/tinc);
};

1.3 Solving an equation¶

Once the equation to be solver has been defined as an instance of the de_system class, the following additional information is necessary to solve it:

initial conditions, \(x(t_i)\) and \(\dot{x}(t_f)\),
the range of integration, from \(t_i\) and \(t_f\),
(optional) set of timepoints at which dense output is required,
(optional) order of WKB approximation to use, order=3,
(optional) relative tolerance, rtol=1e-4,
(optional) absolute tolerance atol=0.0,
(optional) initial step h_0=1,
(optional) output file name full_output="",

Note the following about the optional arguments:

rtol, atol are tolerances on the local error. The global error in the solution is not guaranteed to stay below these values, but the error per step is. In the RK regime (not oscillatory solution), the global error will rise above the tolerance limits, but in the WKB regime, the global error usually stagnates.
The initial step should be thought of as an initial estimate of what the first stepsize should be. The solver will determine the largest possible step within the given tolerance limit, and change h_0 if necessary.
The full output of solve() will be written to the filename contained in full_output, if specified.

Here’s an example to illustrate usage of all of the above variables:

#include "solver.hpp" // de_system, Solution defined in here

// Define the system
de_system sys(...) // For args see previous examples

// Necessary parameters:
// initial conditions
std::complex<double> x0=std::complex<double>(1.0,1.0), dx0=0.0;
// range of integration
double ti=1.0, tf=100.0;

// Optional parameters:
// dense output will be required at the following points:
int n = 1000;
std::vector t_eval(n);
for(int i=0; i<n; i++){
    t_eval[i] = i/10.0;
}
// order of WKB approximation to use
int order=2;
// tolerances
double rtol=2e-4, atol=0.0;
// initial step
double h0 = 0.5;
// write the solution to a file
std::string outfile="output.txt";

Solution solution(sys, x0, dx0, ti, tf, t_eval.data(), order, rtol, atol, h0, outfile);
// Solve the equation:
solution.solve()

Here, we’ve also called the solve() method of the Solution class, to carry out the integration. Now all information about the solution is in solution (and written to output.txt).

1.4 Using the solution¶

Let’s break down what solution contains (what Solution.solve() returns). An instance of a Solution object is returned with the following attributes:

times [std::list of double]: timepoints at which the solution was determined. These are not supplied by the user, rather they are internal steps that the solver has takes. The list starts with \(t_i\) and ends with \(t_f\), these points are always guaranteed to be included.
sol [std::list of std::complex<double>]: the solution at the timepoints specified in times.
dsol [std::list of std::complex<double>]: first derivative of the solution at timepoints specified in times.
wkbs [std::list of int/bool]: types of steps takes at each timepoint in times. 1 if the step was WKB, 0 if it was RK.
ssteps [int]: total number of accepted steps.
totsteps [int]: total number of attempted steps (accepted + rejected).
wkbsteps [int]: total number of successful WKB steps.
x_eval [std::list of std::complex<double>]: dense output, i.e. the solution evaluated at the points specified in the t_eval optional argument
dx_eval [std::list of std::complex<double>]: dense output of the derivative of the solution, evaluted at the points specified in t_eval optional argument.