2
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ 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. $\endgroup$
    – Richard
    Feb 21, 2022 at 22:59
  • 2
    $\begingroup$ 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. $\endgroup$
    – Laurent90
    Feb 22, 2022 at 15:04

1 Answer 1

2
$\begingroup$

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
$\endgroup$
1
  • $\begingroup$ This is exactly what I had in mind, thanks! $\endgroup$ Mar 1, 2022 at 6:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.