boost odeint gives very different values from Python3 scipy-CodePudding

I am trying to integrate a very simple ODE using boost odeint. For some cases, the values are the same (or very similar) to Python's scipy odeint function. But for other initial conditions, the values are vastly different.

The function is: d(uhat) / dt = - alpha^2 * kappa^2 * uhat where alpha is 1.0, and kappa is a constant depending on the case (see values below).

I have tried several different ODE solvers from boost, and none seem to work.

In the code below, the first case gives nearly identical results, the 2nd case is kind of trivial (but reassuring), and the 3rd case gives erroneous answers in the C version.

Here is the C version:

#include <boost/numeric/odeint.hpp>
#include <cstdlib>
#include <iostream>

// The rhs:  d_uhat_dt = - alpha^2 * kappa^2 * uhat
class Eq {
public:
  Eq(double alpha, double kappa)
    : m_constant(-1.0 * alpha * alpha * kappa * kappa) {}
  
  void operator()(double uhat, double& d_uhat_dt, const double t) const
  {
    d_uhat_dt = m_constant * uhat;
  }

private:
  double m_constant;
};

void integrate(double kappa, double initValue)
{
  const unsigned numTimeIncrements = 100;
  const double dt = 0.1;
  const double alpha = 1.0;

  double uhat = initValue;      //Init condition
  std::vector<double> uhats;    //Results vector

  Eq rhs(alpha, kappa);         //The RHS of the ODE
  boost::numeric::odeint::runge_kutta_dopri5<double> stepper;
  
  for(unsigned step = 0; step < numTimeIncrements;   step) {
    uhats.push_back(uhat);
    stepper.do_step(rhs, uhat, step*dt, dt);
  }

  std::cout << "kappa = " << kappa << ", initial value = " << initValue << std::endl;
  for(auto val : uhats)
    std::cout << val << std::endl;
  std::cout << "---" << std::endl << std::endl;
}

int main() {

  const double kappa1 = 0.062831853071796;
  const double initValue1 = -187.097241230045967;
  integrate(kappa1, initValue1);

  const double kappa2 = 28.274333882308138;
  const double initValue2 = 0.000000000000;
  integrate(kappa2, initValue2);

  const double kappa3 = 28.337165735379934;
  const double initValue3 = -0.091204068895190;
  integrate(kappa3, initValue3);
  
  return EXIT_SUCCESS;
}

and the corresponding Python3 version:

enter code here
#!/usr/bin/env python3

import numpy as np
from scipy.integrate import odeint

def Eq(uhat, t, kappa, a):
    d_uhat = -a**2 * kappa**2 * uhat
    return d_uhat

def integrate(kappa, initValue):
    dt = 0.1
    t = np.arange(0,10,dt)
    a = 1.0
    print("kappa = "   str(kappa))
    print("initValue = "   str(initValue))
    uhats = odeint(Eq, initValue, t, args=(kappa,a))
    print(uhats)
    print("---")
    print()

kappa1 = 0.062831853071796
initValue1 = -187.097241230045967
integrate(kappa1, initValue1)

kappa2 = 28.274333882308138
initValue2 = 0.000000000000
integrate(kappa2, initValue2)

kappa3 = 28.337165735379934
initValue3 = -0.091204068895190
integrate(kappa3, initValue3)

CodePudding user response：

This has all the hallmarks of precision issues.

Simply replacing double with long double gives:

Live On Compiler Explorer

#include <boost/numeric/odeint.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <fmt/ranges.h>
#include <fmt/ostream.h>
#include <iostream>

using Value = long double;
//using Value = boost::multiprecision::number<
    //boost::multiprecision::backends::cpp_bin_float<100>,
    //boost::multiprecision::et_off>;

// The rhs:  d_uhat_dt = - alpha^2 * kappa^2 * uhat
class Eq {
  public:
    Eq(Value alpha, Value kappa)
        : m_constant(-1.0 * alpha * alpha * kappa * kappa)
    {
    }

    void operator()(Value uhat, Value& d_uhat_dt, const Value) const
    {
        d_uhat_dt = m_constant * uhat;
    }

  private:
    Value m_constant;
};

void integrate(Value const kappa, Value const initValue)
{
    const unsigned numTimeIncrements = 100;
    const Value    dt                = 0.1;
    const Value    alpha             = 1.0;

    Value              uhat = initValue; // Init condition
    std::vector<Value> uhats;            // Results vector

    Eq rhs(alpha, kappa); // The RHS of the ODE
    boost::numeric::odeint::runge_kutta_dopri5<Value> stepper;

    for (unsigned step = 0; step < numTimeIncrements;   step) {
        uhats.push_back(uhat);
        auto&& stepdt = step * dt;
        stepper.do_step(rhs, uhat, stepdt, dt);
    }

    fmt::print("kappa = {}, initial value = {}\n{}\n---\n", kappa, initValue,
               uhats);
}

int main() {
    for (auto [kappa, initValue] :
         {
             std::pair //
              { 0.062831853071796 , -187.097241230045967 },
              {28.274333882308138 ,    0.000000000000    },
              {28.337165735379934 ,   -0.091204068895190 },
         }) //
    {
        integrate(kappa, initValue);
    }
}

Prints

kappa = 0.0628319, initial value = -187.097
[-187.097, -187.023, -186.95, -186.876, -186.802, -186.728, -186.655, -186.581, -186.507, -186.434, -186.36, -186.287, -186.213, -186.139, -186.066, -185.993, -185.919, -185.846, -185.772, -185.699, -185.626, -185.553, -185.479, -185.406, -185.333, -185.26, -185.187, -185.114, -185.04, -184.967, -184.894, -184.821, -184.748, -184.676, -184.603, -184.53, -184.457, -184.384, -184.311, -184.239, -184.166, -184.093, -184.021, -183.948, -183.875, -183.803, -183.73, -183.658, -183.585, -183.513, -183.44, -183.368, -183.296, -183.223, -183.151, -183.079, -183.006, -182.934, -182.862, -182.79, -182.718, -182.645, -182.573, -182.501, -182.429, -182.357, -182.285, -182.213, -182.141, -182.069, -181.998, -181.926, -181.854, -181.782, -181.71, -181.639, -181.567, -181.495, -181.424, -181.352, -181.281, -181.209, -181.137, -181.066, -180.994, -180.923, -180.852, -180.78, -180.709, -180.638, -180.566, -180.495, -180.424, -180.353, -180.281, -180.21, -180.139, -180.068, -179.997, -179.926]
---
kappa = 28.2743, initial value = 0
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
---
kappa = 28.3372, initial value = -0.0912041
[-0.0912041, -38364100, -1.61375e 16, -6.78809e 24, -2.85534e 33, -1.20107e 42, -5.0522e 50, -2.12516e 59, -8.93928e 67, -3.76022e 76, -1.5817e 85, -6.65327e 93, -2.79864e 102, -1.17722e 111, -4.95186e 119, -2.08295e 128, -8.76174e 136, -3.68554e 145, -1.55029e 154, -6.52114e 162, -2.74306e 171, -1.15384e 180, -4.85352e 188, -2.04159e 197, -8.58774e 205, -3.61235e 214, -1.5195e 223, -6.39163e 231, -2.68858e 240, -1.13092e 249, -4.75713e 257, -2.00104e 266, -8.41718e 274, -3.54061e 283, -1.48932e 292, -6.26469e 300, -2.63518e 309, -1.10846e 318, -4.66265e 326, -1.9613e 335, -8.25002e 343, -3.47029e 352, -1.45975e 361, -6.14028e 369, -2.58285e 378, -1.08645e 387, -4.57005e 395, -1.92235e 404, -8.08618e 412, -3.40137e 421, -1.43075e 430, -6.01833e 438, -2.53155e 447, -1.06487e 456, -4.47929e 464, -1.88417e 473, -7.92559e 481, -3.33382e 490, -1.40234e 499, -5.89881e 507, -2.48128e 516, -1.04373e 525, -4.39033e 533, -1.84675e 542, -7.76818e 550, -3.26761e 559, -1.37449e 568, -5.78166e 576, -2.432e 585, -1.023e 594, -4.30314e 602, -1.81008e 611, -7.61391e 619, -3.20272e 628, -1.34719e 637, -5.66684e 645, -2.3837e 654, -1.00268e 663, -4.21768e 671, -1.77413e 680, -7.4627e 688, -3.13911e 697, -1.32044e 706, -5.55429e 714, -2.33636e 723, -9.82768e 731, -4.13392e 740, -1.73889e 749, -7.31449e 757, -3.07677e 766, -1.29421e 775, -5.44399e 783, -2.28996e 792, -9.6325e 800, -4.05182e 809, -1.70436e 818, -7.16922e 826, -3.01567e 835, -1.26851e 844, -5.33587e 852]
---

As you can see, I tried some simple takes to get it to use Boost Multiprecision, but that didn't immediately work and may require someone with actual understanding of the maths / IDEINT.

CodePudding user response：

With boost::odeint you should use the integrate... interface functions.

What happens in your code using do_step is that you use the dopri5 method as a fixed-step method with your provided step size. In connection with the large coefficient L=kappa^2 of about 800, the product L*dt=80 is far outside the stability region of the method, the boundary is between values of 2 to 3.5. Divergence or oscillating divergence is the expected result.

What should happen and what is implemented in the scipy odeint, ode-dopri5 or solve_ivp-RK45 methods is a method with adaptive step size. Internally the optimal step size for the error tolerances is chosen, and in each internal step an interpolation polynomial is determined. The output or observed values are computed with this interpolator, also called "dense output" if the interpolation function is returned from the integrator. One result of the error control is that the step size is always inside the stability region, for stiff problems with medium error tolerances on or close to the boundary of this region.