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

This is a combination of these two previous questions:

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.

• I'd recommend looking for something in Julia's DifferentialEquations.jl and the surrounding ecosystem. It seems more sophisticated and active than Python or C++'s ecosystem for this purpose. Feb 21, 2022 at 22:59
• I guess if the OP is using C++, changing to another language is not an advisable solution for such a simple problem. Using a dense output is indeed the way to go. Alternatively, you can also split the integration in multiple intervals (one interval for each log-space value). Look into odeint's steppers, it seems to be what you are looking for, hwoever I have no experience myself with this C++ library. Feb 22, 2022 at 15:04

## 1 Answer

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

• This is exactly what I had in mind, thanks! Mar 1, 2022 at 6:21