在基类和派生类中分离 ODE 和 ODE 求解器

如何解决在基类和派生类中分离 ODE 和 ODE 求解器

我觉得这个问题有点长，所以我认为最好先考虑一下它的简化版本：

有两个类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 &,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,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},gamma{gamma},state{state}
    { }
    //---------------------------------------------------//
    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 x0{0.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&,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>&,double)’}
   93 |                  &HarmonicOscillator::dampedOscillator,|                  ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      |                  |
      |                  dim1 (HarmonicOscillator::*)(const dim1&,double) {aka std::vector<double> (HarmonicOscillator::*)(const std::vector<double>&,double)}
example.cpp:29:31: note:   initializing argument 2 of ‘void ODE_Solver::integrate(Integrator,derivative,int,double,double)’
   29 |                    derivative ode_system,|                    ~~~~~~~~~~~^~~~~~~~~~

如果我在同一类 ode 求解器中定义了 ode，它会起作用 link to github。 那么你的想法是什么？

解决方法

您的代码有几个问题。首先，当您想将任意函数作为参数传递给其他函数时，请考虑使用 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。