boost odeint 给出了与 Python3 scipy 非常不同的值
Posted
技术标签:
【中文标题】boost odeint 给出了与 Python3 scipy 非常不同的值【英文标题】：boost odeint gives very different values from Python3 scipy 【发布时间】：2021-12-01 22:14:36 【问题描述】：我正在尝试使用 boost odeint 集成一个非常简单的 ODE。 在某些情况下，这些值与 Python 的 scipy odeint 函数相同（或非常相似）。 但对于其他初始条件，值则大不相同。
函数为：d(uhat) / dt = - alpha^2 * kappa^2 * uhat 其中 alpha 为 1.0，kappa 是一个常数，具体取决于具体情况（请参见下面的值）。
我从 boost 中尝试了几种不同的 ODE 求解器，但似乎都没有。
更新：下面的代码现在可以工作了。
在下面的代码中，第一种情况给出了几乎相同的结果，第二种情况有点微不足道（但令人放心），第三种情况在 C++ 版本中给出了错误的答案。
这里是 C++ 版本：
#include <boost/numeric/odeint.hpp>
#include <cstdlib>
#include <iostream>
typedef boost::numeric::odeint::runge_kutta_dopri5<double> Stepper_Type;
struct ResultsObserver
std::ostream& m_out;
ResultsObserver( std::ostream &out ) : m_out( out )
void operator()(const State_Type& x , double t ) const
m_out << t << " : " << x << std::endl;
;
// 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
//This is what I was doing that did not work
//
//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);
//
//This works
integrate_const(
boost::numeric::odeint::make_dense_output<Stepper_Type>( 1E-12, 1E-6 ),
rhs, uhat, startTime, endTime, dt, ResultsObserver(std::cout)
);
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;
以及对应的Python3版本：
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)
【问题讨论】：
结果有什么不同，它们是完全不同的值还是只是在一些数字之后不同？您是否比较了默认误差容限是否大致相同？最好明确设置。 第一种情况仅在精度的第 8 位或第 9 位不同。第二种情况是全 0，在双精度范围内，所以它们是相等的。但是，第 3 种情况相差 100 多个数量级。它最终射向 -1e300，然后进入 NAN，而 Python 版本进入 / 接近 0，即 -1e-16。 对我上面的评论稍作修正：它最终开始振荡，增长到 + 或 -1e300，然后进入 nan，而 Python 版本平稳地下降到 0，结束于 4e-26 左右上一学期。 【参考方案1】：这具有精度问题的所有特征。
只需将double 替换为long double 即可：
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);
打印
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]
---
如您所见，我尝试了一些简单的方法来让它使用 Boost Multiprecision，但这并没有立即奏效，可能需要真正了解数学/IDEINT 的人。
【讨论】：
【参考方案2】：使用 boost::odeint 你应该使用integrate... 接口函数。
在使用do_step 的代码中发生的情况是，您将dopri5 方法用作具有您提供的步长的固定步长方法。与大约 800 的大系数 L=kappa^2 相关，乘积 L*dt=80 远远超出方法的稳定区域，边界在 2 到 3.5 之间。背离或振荡背离是预期的结果。
在 scipy odeint、ode-dopri5 或 solve_ivp-RK45 方法中应该发生什么以及实现什么是具有自适应步长的方法。在内部选择误差容限的最佳步长，并在每个内部步长中确定一个插值多项式。使用此插值器计算输出或观察值，如果插值函数从积分器返回，则也称为“密集输出”。误差控制的一个结果是步长总是在稳定区域内，对于在该区域边界上或接近该区域边界具有中等误差容限的刚性问题。
【讨论】：
是的，这是正确的。我更新了代码以显示它正常工作。谢谢你的回答。以上是关于boost odeint 给出了与 Python3 scipy 非常不同的值的主要内容，如果未能解决你的问题，请参考以下文章
使用Armadillo和boost :: numeric :: odeint进行模板实例化
如何将向量传递给基于推力的 odeint 观察者的构造函数，以便可以在函子中读取它