c++ - C++ 中 ODE 类的设计模式

Question

我一直在 C++ 中使用不同的 ODE 求解器，并且我熟悉 ODE 的标准和更高级数值方法背后的理论。我想了解的是 ODE 类的“设计模式”是什么。例如，查看 this question
我们注意到 r.h.s 的定义由 void 组成。引用 z 的函数和 dzdt如下:

void ode( const state_type &z , state_type &dzdt , double t ) { 
    dzdt[0] = z[1]; 
    dzdt[1] = -1 * z[0] * w * w; 
} 

然后集成主要进行

int main() { ... 

    integrate( ode , z , t , 1000 , 0.1 , write_ode ); 
    return 0;

}

当然，这样的库确实是硬编码的，但我只想掌握“步进器”背后的一般概念，比如显式欧拉方法。

由于 r.h.s 的定义类似于 void ode (...) ，我想在步进部分有一个调用 void ode(...)允许更新 dzdt .它可以实现如下(使用 std::vector 类)

void do_step(std::vector<double>& z, double tn, double h){
    //tn current time, h time step
    std::vector<double> dzdt(2); 
    ode(tn,y,dzdt);
    z[0] += h*dzdt[0];
    z[1] += h*dzdt[1];
}

IE。:

我为新的解 vector 分配空间

计算新状态

更新 z[0] , z[1]使用欧拉的方案，或其他。

可以用 do_step(y,tn,h); 调用
总而言之，我的问题是:给定 r.h.s. 的定义，这是定义数值方法步骤的好方法吗？任何具有此类问题的一些设计技术的引用/书籍都将受到高度赞赏

编辑
我想了解我是否正确理解了 Lutz 答案的第一段:

The first idea is that the derivatives vectors for the stages of the RK method are made components of the solver class. This prevents frequent memory allocation and de-allocation/garbage collection during the run of the integrator.

为了掌握这个想法，我为谐振子 x''=-(k/m)x 写了一个经典的 RK4，并检查了积分的好坏。我不想进行任何代码调试，但我真的不知道这是否是 Lutz 的意思，因为该函数在这里定义为私有(private)成员，这当然不是一个好点。我想使用 struct用于定义 rhs。在类里面，但我不明白以哪种方式。

#include <iostream>
#include <cmath>
#include <vector>


std::vector<double> operator+(const std::vector<double>& a, const std::vector<double>& b){
    std::vector<double> ret(a.size());
    for(unsigned int i=0;i<a.size();++i){
        ret[i] = a[i]+b[i];
    }
    return ret;
}


constexpr double k = 3.0;
constexpr double m = 2.0;
constexpr double km  = k/m;

class Rk{
    
private:
    std::vector<double> f(const double t, const std::vector<double>& y){
        std::vector<double> state(2);
        state[0] = y[1];
        state[1] = -(k/m)*y[0] + t;
        return state;
    }
    
    std::vector<double> y0;
    const double T;
    double dt;
    
public:
    Rk( std::vector<double> _y0,const double _T,double _dt) : y0{_y0}, T{_T}, dt{_dt}{}
    ~Rk()=default;
    
    
    
    std::vector<double> mvec(const std::vector<double>& v, const double c) {
        //implements m*vec
        const auto size = v.size();
        std::vector<double> res(size);
        for (std::size_t i = 0; i < size; ++i){
            res[i] = c*v[i];
        }
        return res;
    }
    
    
    
    
    void step(std::vector<double>& state, const double t){
        //performs a Rk4 step from time t to t+dt
        //state: current state y_1(tn),y_2(tn),...
        const double dth = 0.5*dt;
        std::vector<double> k1 = f(t,state);
        std::vector<double> k2 = f(t+dth,state+mvec(k1,0.5*dt));
        std::vector<double> k3 = f(t+dth,state+mvec(k2,0.5*dt));
        std::vector<double> k4 = f(t+dt,state+mvec(k3,dt));
        state = state + mvec ((k1+mvec(k2,2)+mvec(k3,2) + k4),dt/6.0);
    }
    
    void integrate(){
        std::vector<double> state(2);
        state = y0;
        double t =  0.0;
        for (unsigned int i=0;i<std::ceil(T/dt);++i){
            step(state,t);
            t+=dt;
            //            double err = std::fabs(std::sqrt(1/km) * std::sin(std::sqrt(km)*(t)) - state[0]);
            double err = std::fabs((1.0/9.0) * (6*t+std::sqrt(6)*std::sin(std::sqrt(km)*t)) - state[0]);
            std::cout << err <<std::endl;
            
        }
    }
    
    
};







int main(){
    const double dt = 0.01;
    std::vector<double> y0;
    y0.push_back(0.0);
    y0.push_back(1.0);
    Rk my_ode{y0,1.0,dt};
    my_ode.integrate();
    
    
    return 0;
    
}

Answer 1

第一个想法是将 RK 方法阶段的导数 vector 作为求解器类的组件。这可以防止在积分器运行期间频繁的内存分配和取消分配/垃圾收集。使用高阶方法可能会更好，看看为什么这很有用，欧拉太琐碎了，可能会给出错误的直觉。
下一个要考虑的想法是使用可变的、适应的步长方法。这意味着您可以在需要时执行内部步骤，并通过插值评估外部使用，通常使用“密集输出”概念。在那里您可以完全隐藏内部步骤，或公开它们和插值函数/对象。这两个想法都可以在 scipy.integrate 中进行研究。求解器，旧的步进器类 ode隐藏内部步骤，步进器类 RK45, Radau,...背后新solve_ivp接口(interface)实现第二个概念。
然后你很快就到了你有一个结构化的状态空间的地步。可以实现这样一种理念，即数据必须以扁平的一维 vector 组合，这是(旧的)标准。或者可以为状态空间类配备必要的 vector 和范数运算。您应该已经在 boost::odeint 模板参数中发现了这一点。
下一点是，状态空间可能有不同的对象、位置与速度等部分，它们的大小尺度差异很大，应该在误差估计(包括绝对误差容限)中单独处理步长 Controller (即，计算每个段的最佳步长，然后取最小值)。
最后一点，求解器只有在具有事件- Action 机制时才会变得普遍有用。事件是状态的某些功能的(定向)零交叉， Action 可以是记录事件，终止事件，或者不太标准，对状态 vector 的修改。