C ++中ODE类的设计模式 [英] Design pattern of an ODE class in C++
问题描述
我一直在C ++中使用不同的ODE求解器，并且我熟悉ODE的标准和更高级的数值方法背后的理论.我想了解的是什么是设计模式"ODE类.例如，查看此问题
I've been using different ODE solvers in C++, and I'm familiar with the theory behind standard and more advanced numerical methods for ODEs. What I'd like to understand is what is the "design pattern" of an ODE class. For instance, looking at this question
我们注意到rhs的定义由 void 函数组成，该函数采用对 z 和 dzdt 的引用，如下所示:
we notice that the definition of the r.h.s consists of a void function which takes references to z and dzdt as follows:
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;
}
当然，这样的库确实是硬编码的，但是我只想掌握步进器"背后的一般思想，比如说明确的欧拉方法.
Of course such library is really hard-coded, but I just want to grasp the general idea behind the "stepper", let's say for explicit Euler's method.
由于rhs的定义类似于 void ode(...)，因此我想在步进器部分调用了 void ode(...)可以更新 dzdt .可以如下实现(使用 std :: vector 类)
Since the r.h.s is defined like in void ode (...), I imagine that in the stepper part there's a call to void ode(...) that allows to update dzdt. It could be implemented as follows (using the std::vector class)
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];
}
即:
- 我为新的求解向量分配空间
- 计算新状态
- 使用欧拉方案或其他方法更新
z [0]，z [1].
可以用 do_step(y，tn，h);
总结一下，我的问题是:给定r.h.s.的定义，这是定义数值方法步骤的好方法吗?任何具有这种问题的某些设计技术的参考书/书都将受到高度赞赏
Summarizing, my question is: given the definition of the r.h.s., is this the good way to define a step of a numerical method? Any reference/book with some design techniques of such a problem is highly appreciated
编辑
我想了解一下我是否正确理解了Lutz的答案的第一段:
I'd like to understand if I got correctly the first paragraph of Lutz's answer:
第一个想法是导数阶的导数向量RK方法是求解器类的组成部分.这样可以防止期间频繁的内存分配和取消分配/垃圾收集积分器的运行.
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的意思，因为在这里将该函数定义为私有成员，这当然不是一个好主意.我想使用 struct 定义rhs.在课堂上，但我不知道用哪种方式.
In order to grasp the idea I wrote a classical RK4 for the harmonic oscillator x''=-(k/m)x, and also checked the goodness of the integration. I don't want any code debugging, but I don't really know if it was Lutz's meant, because the function is defined as a private member here which is of course not a good point. I'd like to use struct for defining the rhs. inside the class, but I don't understand in which way.
#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;
}
推荐答案
第一个想法是，将RK方法各阶段的导数向量作为Solver类的组成部分.这样可以防止在集成器运行期间频繁进行内存分配和取消分配/垃圾收集.最好使用高阶方法来查看为什么这样做很有用，因为欧拉太琐碎，可能会给出错误的直觉.
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. It would perhaps be better to use a higher-order method to see why this is useful, Euler is too trivial, may give a false intuition.
要考虑的下一个想法是使用可变的，经过调整的步长方法.这意味着您具有在需要时执行的内部步骤，并通过插值(通常使用密集输出")进行评估以供外部使用.概念.在这里，您可以完全隐藏内部步骤，也可以公开它们以及插值函数/对象.可以在 scipy.integrate 求解器中研究这两种思想，旧的步进器类 ode 隐藏了内部步骤，步进器类 RK45，Radau，... <新的 solve_ivp 接口后面的/code>实现了第二个概念.
The next idea to contemplate is the use of variable, adapted step size methods. This means that you have internal steps that are performed when needed, and evaluation for external use via interpolation, usually using the "dense output" concept. There you could hide the internal steps completely, or expose them and the interpolation function/object. Both ideas can be studied in the scipy.integrate solvers, the old stepper class ode hides the internal steps, the stepper classes RK45, Radau,... behind the new solve_ivp interface implement the second concept.
然后，您迅速到达拥有结构化状态空间的地步.可以实现这样一种哲学，即数据必须以平坦的一维向量(这是(旧)标准)进行组合.或者可以为状态空间类配备必要的向量和范数运算.您应该已经在boost :: odeint模板参数中找到了.
Then you come rapidly to the point where you have a structured state space. One could implement the philosophy that the data has to be assembled in a flat one-dimensional vector, which is (the old) standard. Or one could equip the state space class with the necessary vector and norm operations. You should have found that in the boost::odeint template parameters.
接下来的一点是，状态空间可能具有段，例如不同的对象，位置与速度等，它们的大小比例差异很大，应在误差估计(包括绝对误差容限)中分别进行处理.)作为步长控制器(也就是说，为每个细分计算最佳步长，然后取最小值).
A next point to that is that the state space may have segments, like different objects, position vs. velocity, etc., that are on vastly different magnitude scales that should be treated separately in error estimates (incl. absolute error tolerances) for the step size controller (that is, compute optimal step sizes for every segment and then take the minimum).
最后一点是，只有具有事件动作机制的求解器才变得普遍有用.事件是状态的某些功能的(定向)零交叉，动作可以是记录事件，在事件处终止或以不太标准的方式对状态向量进行修改.
As a last point here, a solver only becomes universally useful if it has an event-action mechanism. Events are (directional) zero-crossings of some functions of the state, actions can be to record the event, to terminate at the event, or less standard, a modification of the state vector.
这篇关于C ++中ODE类的设计模式的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
