在基类和派生类中分离 ODE 和 ODE 求解器
Posted
技术标签:
【中文标题】在基类和派生类中分离 ODE 和 ODE 求解器【英文标题】：Seperating ODE and ODE solver in base and derived classes 【发布时间】：2021-06-13 13:33:45 【问题描述】：我觉得这个问题有点长，所以我认为最好先考虑一下它的简化版本：
有两个类A和B。B继承自A。B中有一个成员函数（add）需要使用A中的成员函数运行。
class A;
typedef int(A::*operation)(int c);
typedef void (A::*myfunction)(operation, int);
class A
public:
int a;
int b;
int do_something(myfunction f, operation op)
(this->*f)(op, 1);
void dummy_name(operation op, int a)
int c = (this->*op)(a);
;
class B : public A
public:
int a, b;
B(int a, int b): a(a), b(b)
int add(int c)
return a+b+c;
;
int main()
B inst_B(2, 5);
inst_B.do_something(&A::dummy_name, &B::add);
simple.cpp:45:41: error: cannot convert ‘int (B::*)(int)’ to ‘operation’ aka ‘int (A::*)(int)’
45 | inst_B.do_something(&A::dummy_name, &B::add);
| ^~~~~~~
| |
| int (B::*)(int)
simple.cpp:17:47: note: initializing argument 2 of ‘void A::do_something(myfunction, operation)’
17 | void do_something(myfunction f, operation op)
| ~~~~~~~~~~^~
为了编写一个简单的 ode 求解器并避免在模型类中使用积分器，对于每个包含常微分方程组的模型，我将求解器和方程分为两类，而模型继承自 ode 求解器。
class HarmonicOscillator: public ODE_Sover
这是一个简化的示例，其中包含一些参数。为了避免传递许多参数和抽象，我更喜欢在一个类中定义 ODE。
我还为导数（dy/dt = f'(y) 的右侧）和积分器（这里只有欧拉积分器）使用了两个函数模板。 这是我想出的：
#include <iostream>
#include <assert.h>
#include <random>
#include <vector>
#include <string>
using std::string;
using std::vector;
class ODE_Solver;
class HarmonicOscillator;
typedef vector<double> dim1;
typedef dim1 (ODE_Solver::*derivative)(const dim1 &, dim1&, const double t);
typedef void (ODE_Solver::*Integrator)(derivative, dim1 &, dim1&, const double t);
class ODE_Solver
public:
ODE_Solver()
double t;
double dt;
dim1 state;
dim1 dydt;
void integrate(Integrator integrator,
derivative ode_system,
const int N,
const double ti,
const double tf,
const double dt)
dim1 dydt(N);
const size_t num_steps = int((tf-ti) / dt);
for (size_t step = 0; step < num_steps; ++step)
double t = step * dt;
(this->*integrator)(ode_system, state, dydt, t);
// print state
void eulerIntegrator(derivative RHS, dim1 &y, dim1 &dydt, const double t)
int n = y.size();
(this->*RHS)(y, dydt, t);
for (int i = 0; i < n; i++)
y[i] += dydt[i] * dt;
;
class HarmonicOscillator: public ODE_Solver
public:
int N;
double dt;
double gamma;
string method;
dim1 state;
// constructor
HarmonicOscillator(int N,
double gamma,
dim1 state
) : N N, gammagamma, statestate
//---------------------------------------------------//
dim1 dampedOscillator(const dim1 &x, dim1&dxdt, const double t)
dxdt[0] = x[1];
dxdt[1] = -x[0] - gamma * x[1];
return dxdt;
;
//-------------------------------------------------------//
int main(int argc, char **argv)
const int N = 2;
const double gamma = 0.05;
const double t_iinit = 0.0;
const double t_final = 10.0;
const double dt = 0.01;
dim1 x00.0, 1.0;
HarmonicOscillator ho(N, gamma, x0);
ho.integrate(&ODE_Solver::eulerIntegrator,
&HarmonicOscillator::dampedOscillator,
N, t_iinit, t_final, dt);
return 0;
我收到以下错误：
example.cpp: In function ‘int main(int, char**)’:
example.cpp:93:18: error: cannot convert ‘dim1 (HarmonicOscillator::*)(const dim1&, dim1&, double)’ aka ‘std::vector<double> (HarmonicOscillator::*)(const std::vector<double>&, std::vector<double>&, double)’ to ‘derivative’ aka ‘std::vector<double> (ODE_Solver::*)(const std::vector<double>&, std::vector<double>&, double)’
93 | &HarmonicOscillator::dampedOscillator,
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
| |
| dim1 (HarmonicOscillator::*)(const dim1&, dim1&, double) aka std::vector<double> (HarmonicOscillator::*)(const std::vector<double>&, std::vector<double>&, double)
example.cpp:29:31: note: initializing argument 2 of ‘void ODE_Solver::integrate(Integrator, derivative, int, double, double, double)’
29 | derivative ode_system,
| ~~~~~~~~~~~^~~~~~~~~~
如果我在同一类 ode 求解器中定义 ode，它可以工作 link to github。 那么你的想法是什么？
【问题讨论】：
你的类程序集缠绕得太紧，减少连接，给求解器和IVP良好分离的类。您可以在 ***.com/questions/49738918/solving-lorenz-equation 中找到这样的设置，对于其他不只是 C 伪装成 C++ 的方法，请参阅我在 ***.com/questions/59231161/…、***.com/questions/34965829/runge-kuttas-method 中的回答 您还可以研究 boost-odeint 或 julia 包中的接口以获取灵感，但在这些扩展功能中，例如步长控制、密集输出、事件-动作机制，...对于固定步长的显式 Runge-Kutta 方法，并非所有这些复杂性都是必需的。 感谢 cmets，我已经查看了 boost.odeint，但我也想尝试这种方法，因为它的经验。 odeint 需要将参数传递给每次调用的导数。我还需要定义一些 odeint 不支持的 DDE 求解器。这只是一个尝试。 【参考方案1】：您的代码存在几个问题。首先，当您想将任意函数作为参数传递给其他函数时，请考虑使用std::function。
其次，应该使用继承来声明"is-a" relationships。由于谐波振荡器不是 ODE 求解器，因此不要为此使用继承。也没有“has-a”关系（所以composition 也不合适），而是求解器作用于给定函数，所以最合适的做法是传递谐波振荡器作为求解器函数的参数。
代码可能看起来的示例：
class HarmonicOscillator
...
public:
...
double operator()(double t)
...
return /* value at time t */;
;
double integrate(std::function<double(double)> func, double start, double end, double dt)
double sum = 0;
for (double t = start; t < end; t += dt)
sum += func(t) * dt;
return sum;
然后你就这样称呼它：
HarmonicOscillator ho(...);
auto result = integrate(ho, t_iinit, t_final, dt);
上面的代码可能不是你想要的，但这是我认为你应该瞄准的代码结构。
如果您希望能够处理不仅采用 double 并返回 double 的函数，而且还希望能够处理任意类型，您可以将 integrate() 设为模板：
template <typename Function, typename T>
auto integrate(Function func, T start, T end, T dt)
decltype(func(start)) sum;
for (T t = start; t < end; t += dt)
sum += func(t) * dt;
return sum;
如果您为支持算术运算的输入和输出值创建适当的类型，这将有效，它不适用于您的dim1。我建议你尝试找到一个实现数学向量类型的库，例如Eigen。
【讨论】：
谢谢。对我来说，两个主要目标是避免为每个模型编写积分器（ODE、DDE、SDE、...求解器）以及易于传递多个参数（这里只有一个：gamma，有时最多 20 个）。我考虑了一下。 @Abolfazl ：但这就是 ODE 函数类的意义所在，它的 () 运算符。您可以将所有参数作为该类的数据成员，通过构造函数或特定的 setter 函数设置它们的值，并在 () 运算符中具有 ODE 函数的标准接口。 /// ODE、DDE和SDE在概念上需要不同的接口，没有必要强制它们变成统一的格式。以上是关于在基类和派生类中分离 ODE 和 ODE 求解器的主要内容，如果未能解决你的问题，请参考以下文章