【发布时间】：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 左右上一学期。
标签： c++ python-3.x boost ode odeint