How to use an adaptive step size in boost::odeint

Question

This is a combination of these two previous questions:

• How to get ODE solution at specified time points?
• Stop integration in odeint with stiff ode

I need to solve the following differential equation for $f_N$,

$$\frac{d}{d t} f_{N}(\mathbf{k}, t)=D(\mathbf{k}, t)\left[f_{N}(\mathbf{k}, t)-f_{N}^{\mathrm{eq}}(\mathbf{k}, t)\right]$$

where $\mathbf k$ is a parameter and $t$ is time. I need the values of $f_N$ at specific times $t=T$; where $T$ is a log-spaced array from $10^{-6}$ to $10^3$. In Python, as suggested in the first link, I would do something like this:

solve_ivp(lambda t, y: ODE__system(t, y), 
                     t_span=(1e-6, 1e3), y0=[0], t_eval=T, vectorized=False, 
                     method='BDF', rtol=1e-12, atol=1e-15)

I need to do the same in C++. For the log-spaced times, I have the following piece of code:

std::vector<double> logspace(double a, double b, int k) {
  /*
  y = linspace(start, stop, num=num, endpoint=endpoint, axis=axis)
  if dtype is None:
      return _nx.power(base, y)
  return _nx.power(base, y).astype(dtype, copy=False)
  */
  const auto exp_scale = (b - a) / (k - 1);
  std::vector<double> logspace;
  logspace.reserve(k);
  for (int i = 0; i < k; i++) {
    logspace.push_back(a + i * exp_scale);
  }
  std::for_each(logspace.begin(), logspace.end(),
                [](double &x) { x = pow(10, x); });
  return logspace;
}

...


std::vector<double> T = logspace(-6, 3, 500);

Regarding the ODE, I thought boost::odeint with a dense output stepper would help -as suggested in the second question- but I don't understand how to implement a non-constant step. Any guidance would be highly appreciated.

Answer 1

Here's an example of using Boost odeint with output at predefined times using the integrate_times integration function. (I hardcoded a vector of times {0.0, 1.0, 2.0, 4.0, ...}, but you could use your logspace function to generate the times.) The differential equation is the logistic differential equation.

#include <cstdio>
#include <vector>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

typedef boost::numeric::ublas::vector<double> state_type;


class logistic
{
    double r;  // Small population growth rate
    double K;  // Carrying capacity

public:

    logistic(double r, double K) : r(r), K(K) {}

    void operator()(const state_type &y, state_type &dydt, const double /* t */)
    {
        dydt[0] = r*y[0]*(1 - y[0]/K);
    }
};


void writer(const state_type &y, const double t)
{
    printf("%7.1f ", t);
    for (auto v : y) {
        printf(" %10.7f", v);
    }
    printf("\n");
}


int main(int argc, char *argv[])
{
    state_type y(1);
    double initial_step = 1e-6;
    vector<double> times = {0.0, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0,
                            64.0, 128.0, 256.0, 512.0};

    // Initial condition
    y[0] = 1e-4;

    bulirsch_stoer_dense_out<state_type> stepper(1e-9, 1e-9, 1.0, 0.0);
    auto sys = logistic(0.08, 3.0);
    integrate_times(stepper, sys, y,
                    times.begin(), times.end(),
                    initial_step, writer);
}

Output:

    0.0   0.0001000
    1.0   0.0001083
    2.0   0.0001174
    4.0   0.0001377
    8.0   0.0001896
   16.0   0.0003596
   32.0   0.0012931
   64.0   0.0166413
  128.0   1.4483309
  256.0   2.9998852
  512.0   3.0000000