我问过OP他们是否能够使用模板，他们回答说"I can use templates, but I wanted it to be simple."首先在幕后设置模板结构还需要做一些工作；但一旦完成；模板的使用确实大大简化了事情。最终，它通常会使代码更加通用，而不是特定于单个任务。
这是一个未实现的类，用于显示 ODE 积分器的通用模式。
#ifndef ODE_GENERALIZED_DEFINITION
#define ODE_GENERALIZED_DEFINITION
// A generalized structure of what an ODE should represent
// This class is not used; but only serves to show the interface
// of what any ODE type integrator - solver needs to contain
template<class Container, class Time = double, class Traits = container_traits<Container> >
class ode_step {
public:
typedef unsigned short order_type;
typedef Time time_type;
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
ode_step() = default;
order_type order_step() const {};
void adjust_size( const container_type& x ) {}
// performs one step
template<class DynamicSystem>
void do_step( DynamicSystem& system, container_type& x, const container_type& dxdt, time_type t, time_type dt ) {}
// performs one step
template<class DynamicSystem>
void do_step( DynamicSystem& system, container_type& x, time_type t, time_type dt ) {}
};
#endif // !ODE_GENERALIZED_DEFINITION
上面不会做任何事情，而仅用于展示 ODE 如何表现或起作用。在我们按照上述模式对任何类型的 ODE 进行建模之前，我们确实需要定义一些东西。
首先我们需要一个container_traits 类。
container_traits.h
#ifndef CONTAINER_TRAITS
#define CONTAINER_TRAITS
template<class Container>
struct container_traits {
typedef Container container_type;
typedef typename container_type::value_type value_type;
typedef typename container_type::iterator iterator;
typedef typename container_type::const_iterator const_iterator;
static void resize( const container_type& x, container_type& dxdt ) {
dxdt.resize( x.size() );
}
static bool same_size( const container_type& x1, const container_type& x2 ) {
return (x1.size() == x2.size());
}
static void adjust_size( const container_type& x1, container_type& x2 ) {
if( !same_size( x1, x2 ) ) resize( x1, x2 );
}
static iterator begin( container_type& x ) {
return x.begin();
}
static const_iterator begin( const container_type& x ) {
return x.begin();
}
static iterator end( container_type& x ) {
return x.end();
}
static const_iterator end( const container_type& x ) {
return x.end();
}
};
#endif // !CONTAINER_TRAITS
上面的类相当简单且直接。
然后我们需要一组迭代器类型函数来迭代这些容器。迭代器函数稍微复杂一些，但它们的预期目的又相当简单。我将它们包装在命名空间中，这样它们就不会与可能存在的任何其他迭代器类或函数发生冲突。
iterator_algebra.h
#ifndef ITERATOR_ALGEBRA
#define ITERATOR_ALGEBRA
#include <cmath>
#include <iostream>
namespace it_algebra { // iterator algebra
// computes y += alpha * x1
template<class InOutIter, class InIter, class T>
void increment( InOutIter first1, InOutIter last1,
InIter first2, T alpha ) {
while( first1 != last1 )
(*first1++) += alpha * (*first2++);
}
// computes y = alpha1 * ( x1 + x2 + alpha2*x3 )
template<class OutIter, class InIter1, class InIter2, class InIter3, class T>
void increment_sum_sum( OutIter first1, OutIter last1,
InIter1 first2, InIter2 first3,
InIter3 first4, T alpha1, T alpha2 ) {
while( first1 != last1 )
(*first1++) += alpha1 *
((*first2++) + (*first3++) + alpha2 * (*first4++));
}
// computes y = alpha1*x1 + alpha2*x2
template<class OutIter, class InIter1, class InIter2, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin ) {
while( y_begin != y_end ) {
(*y_begin++) = alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++);
}
}
// computes y = x1 + alpha2*x2 + alpha3*x3
template<class OutIter, class InIter1, class InIter2, class InIter3, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
template<class OutIter, class InIter1, class InIter2,
class InIter3, class InIter4, class InIter5, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin,
T alpha5, InIter5 x5_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++) +
alpha5 * (*x5_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6
template<class OutIter, class InIter1, class InIter2,
class InIter3, class InIter4, class InIter5, class InIter6, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin,
T alpha5, InIter5 x5_begin,
T alpha6, InIter6 x6_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++) +
alpha5 * (*x5_begin++) +
alpha6 * (*x6_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7
template<class OutIter, class InIter1, class InIter2, class InIter3,
class InIter4, class InIter5, class InIter6, class InIter7, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin,
T alpha5, InIter5 x5_begin,
T alpha6, InIter6 x6_begin,
T alpha7, InIter7 x7_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++) +
alpha5 * (*x5_begin++) +
alpha6 * (*x6_begin++) +
alpha7 * (*x7_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7 + alpha8*x8
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4,
class InIter5, class InIter6, class InIter7, class InIter8, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin,
T alpha5, InIter5 x5_begin,
T alpha6, InIter6 x6_begin,
T alpha7, InIter7 x7_begin,
T alpha8, InIter8 x8_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++) +
alpha5 * (*x5_begin++) +
alpha6 * (*x6_begin++) +
alpha7 * (*x7_begin++) +
alpha8 * (*x8_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7 + alpha8*x8 + alpha9*x9
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4,
class InIter5, class InIter6, class InIter7, class InIter8, class InIter9, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin,
T alpha5, InIter5 x5_begin,
T alpha6, InIter6 x6_begin,
T alpha7, InIter7 x7_begin,
T alpha8, InIter8 x8_begin,
T alpha9, InIter9 x9_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++) +
alpha5 * (*x5_begin++) +
alpha6 * (*x6_begin++) +
alpha7 * (*x7_begin++) +
alpha8 * (*x8_begin++) +
alpha9 * (*x9_begin++);
}
// computes y = x1 + alpha2*x2 + alpha3*x3 + alpha4*x4 + alpha5*x5
// + alpha6*x6 + alpha7*x7 + alpha8*x8 + alpha9*x9 + alpha10*x10
template<class OutIter, class InIter1, class InIter2, class InIter3, class InIter4, class InIter5,
class InIter6, class InIter7, class InIter8, class InIter9, class InIter10, class T>
inline void scale_sum( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin,
T alpha5, InIter5 x5_begin,
T alpha6, InIter6 x6_begin,
T alpha7, InIter7 x7_begin,
T alpha8, InIter8 x8_begin,
T alpha9, InIter9 x9_begin,
T alpha10, InIter10 x10_begin ) {
while( y_begin != y_end )
(*y_begin++) =
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++) +
alpha5 * (*x5_begin++) +
alpha6 * (*x6_begin++) +
alpha7 * (*x7_begin++) +
alpha8 * (*x8_begin++) +
alpha9 * (*x9_begin++) +
alpha10 * (*x10_begin++);
}
// generic version for n values
template<class OutIter, class InIter, class InIterIter, class FactorIter, class T>
inline void scale_sum_generic( OutIter y_begin, OutIter y_end,
FactorIter alpha_begin, FactorIter alpha_end,
T beta, InIter x_begin, InIterIter x_iter_begin ) {
FactorIter alpha_iter;
InIterIter x_iter_iter;
while( y_begin != y_end ) {
x_iter_iter = x_iter_begin;
alpha_iter = alpha_begin;
*y_begin = *x_begin++;
//std::clog<<(*y_begin);
while( alpha_iter != alpha_end ) {
//std::clog<< " + " <<beta<<" * "<<*alpha_iter<<"*"<<(*(*(x_iter_iter)));
(*y_begin) += beta * (*alpha_iter++) * (*(*x_iter_iter++)++);
}
//std::clog<<" = "<<*y_begin<<std::endl;
y_begin++;
}
//std::clog<<std::endl;
}
// computes y += alpha1*x1 + alpha2*x2 + alpha3*x3 + alpha4*x4
template<class OutIter, class InIter1, class InIter2,
class InIter3, class InIter4, class T>
inline void scale_sum_inplace( OutIter y_begin, OutIter y_end,
T alpha1, InIter1 x1_begin,
T alpha2, InIter2 x2_begin,
T alpha3, InIter3 x3_begin,
T alpha4, InIter4 x4_begin ) {
while( y_begin != y_end )
(*y_begin++) +=
alpha1 * (*x1_begin++) +
alpha2 * (*x2_begin++) +
alpha3 * (*x3_begin++) +
alpha4 * (*x4_begin++);
}
// calculates tmp = y, y = x1 + alpha*x2, x1 = tmp
template<class OutIter, class InIter, class T>
inline void scale_sum_swap( OutIter y_begin, OutIter y_end,
OutIter x1_begin, T alpha, InIter x2_begin ) {
T swap;
while( y_begin != y_end ) {
swap = (*x1_begin) + alpha * (*x2_begin++);
*x1_begin++ = *y_begin;
*y_begin++ = swap;
}
}
// computes y = x1 + alpha2 * x2 ; x2 += x3
template<class OutIter, class InIter1,
class InOutIter, class InIter2, class T>
void assign_sum_increment( OutIter first1, OutIter last1, InIter1 first2,
InOutIter first3, InIter2 first4, T alpha ) {
while( first1 != last1 ) {
(*first1++) = (*first2++) + alpha * (*first3);
(*first3++) += (*first4++);
}
}
template<class OutIter, class InIter1, class InIter2, class T >
void weighted_scale( OutIter y_begin, OutIter y_end, InIter1 x1_begin, InIter2 x2_begin,
T eps_abs, T eps_rel, T a_x, T a_dxdt ) {
using std::abs;
while( y_begin != y_end ) {
*y_begin++ = eps_abs +
eps_rel * (a_x * abs( *x1_begin++ ) +
a_dxdt * abs( *x2_begin++ ));
}
}
template<class InIter1, class InIter2, class T >
T max_ratio( InIter1 x1_begin, InIter1 x1_end, InIter2 x2_begin, T initial_max ) {
using std::abs;
while( x1_begin != x1_end ) {
initial_max = std::max( static_cast<T>(abs( *x1_begin++ ) / abs( *x2_begin++ )), initial_max );
}
return initial_max;
}
} // namespace it_algebra
#endif // !ITERATOR_ALGEBRA
现在我们已经有了所需的结构和函数，我们可以使用上面的模式来实现不同类型的 ODE 积分器。我将展示两个最常见的：euler 和 rk4。只要遵循模式，就可以轻松实现 midpoint、rk5 或任何其他方法。
euler.h
#ifndef EULER_H
#define EULER_H
#include "iterator_algebra.h"
#include "container_traits.h"
template<class Container, class Time = double, class Traits = container_traits<Container> >
class euler_stepper {
public:
typedef unsigned short order_type;
typedef Time time_type;
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
private:
container_type m_dxdt;
public:
euler_stepper() = default;
euler_stepper( const container_type& x ) {
adjust_size( x );
}
order_type order_step() const {
return 1;
}
void adjust_size( const container_type& x ) {
traits_type::adjust_size( x, m_dxdt );
}
// performs one step with the knowledge of dxdt(t)
template<class DynamicSystem>
void do_step( DynamicSystem& system, container_type& x, const container_type& dxdt, time_type t, time_type dt ) {
it_algebra::increment( traits_type::begin( x ),
traits_type::end( x ),
traits_type::begin( dxdt ),
dt );
}
// performs one step
template<class DynamicSystem>
void do_step( DynamicSystem& system, container_type& x, time_type t, time_type dt ) {
system( x, m_dxdt, t );
do_step( system, x, m_dxdt, t, dt );
}
};
#endif // EULER_H
rk4.h
#ifndef RK4_H
#define RK4_H
#include "iterator_algebra.h"
#include "container_traits.h"
template<class Container, class Time = double, class Traits = container_traits<Container>>
class rk4_stepper {
public:
typedef unsigned short order_type;
typedef Time time_type;
typedef Traits traits_type;
typedef typename traits_type::container_type container_type;
typedef typename traits_type::value_type value_type;
// typedef typename traits_type::iterator iterator;
// typedef typename traits_type::const_iterator const_iterator;
private:
container_type m_dxdt;
container_type m_dxt;
container_type m_dxm;
container_type m_dxh;
container_type m_xt;
public:
rk4_stepper() = default;
rk4_stepper( const container_type& x ) { adjust_size( x ); }
order_type order_step() const { return 4; }
void adjust_size( const container_type& x ) {
traits_type::adjust_size( x, m_dxdt );
traits_type::adjust_size( x, m_dxt );
traits_type::adjust_size( x, m_dxm );
traits_type::adjust_size( x, m_xt );
traits_type::adjust_size( x, m_dxh );
}
template<class DynamicSystem>
void do_step( DynamicSystem& system, container_type& x, const container_type& dxdt, time_type t, time_type dt ) {
using namespace it_algebra;
const time_type val1 = static_cast<time_type>(1.0);
time_type dh = static_cast<time_type>(0.5) * dt;
time_type th = t + dh;
// dt * dxdt = k1
// m_xt = x + dh*dxdt
scale_sum( traits_type::begin( m_xt ), traits_type::end( m_xt ),
val1, traits_type::begin( x ), dh, traits_type::begin( dxdt ) );
// dt * m_dxt = k2
system( m_xt, m_dxt, th );
// m_xt = x + dh*m_dxt
scale_sum( traits_type::begin( m_xt ), traits_type::end( m_xt ),
val1, traits_type::begin( x ), dh, traits_type::begin( m_dxt ) );
// dt * m_dxm = k3
system( m_xt, m_dxm, th );
// m_xt = x + dt*m_dxm
scale_sum( traits_type::begin( m_xt ), traits_type::end( m_xt ),
val1, traits_type::begin( x ), dt, traits_type::begin( m_dxm ) );
// dt * m_dxh = k4
system( m_xt, m_dxh, t + dt );
// x += dt / 6 * (m_dxdt + m_dxt + val2*m_dxm)
time_type dt6 = dt / static_cast<time_type>(6.0);
time_type dt3 = dt / static_cast<time_type>(3.0);
scale_sum_inplace( traits_type::begin( x ), traits_type::end( x ),
dt6, traits_type::begin( dxdt ),
dt3, traits_type::begin( m_dxt ),
dt3, traits_type::begin( m_dxm ),
dt6, traits_type::begin( m_dxh ) );
}
template<class DynamicSystem>
void do_step( DynamicSystem& system, container_type& x, time_type t, time_type dt ) {
system( x, m_dxdt, t );
do_step( system, x, m_dxdt, t, dt );
}
};
#endif // !RK4_H
现在我们已经拥有了我们所需要的一切；剩下要做的就是使用它们。
OP 提到解决洛伦兹问题，所以我将用它来演示如何使用 ODE：
main.cpp
#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include "euler.h"
#include "rk4.h"
#define tab "\t"
typedef std::vector<double> state_type;
const double sigma = 10.0;
const double R = 28.0;
const double b = 8.0 / 3.0;
void lorenz( state_type& x, state_type& dxdt, double t ) {
dxdt[0] = sigma * (x[1] - x[0]);
dxdt[1] = R * x[0] - x[1] - x[0] * x[2];
dxdt[2] = x[0] * x[1] - b * x[2];
}
int main() {
const double dt = 0.01;
state_type x1( 3 );
x1[0] = 1.0;
x1[1] = 0.0;
x1[2] = 0.0;
state_type x2( 3 );
x2[0] = 1.0;
x2[1] = 0.0;
x2[2] = 0.0;
euler_stepper<state_type> stepper_euler;
stepper_euler.adjust_size( x1 );
rk4_stepper<state_type> stepper_rk4;
stepper_rk4.adjust_size( x2 );
std::fstream file( "compare.txt", std::ios::out );
file << tab << "Euler Stepper to Solve Lorenz" << tab << tab << tab << tab << "RK4 Stepper to Solve Lorenz\n"
<< tab << "========================" << tab << tab << tab << "=======================\n";
double t1 = 0.0;
double t2 = 0.0;
for( size_t oi = 0; oi < 10000; ++oi, t1 += dt, t2 += dt ) {
stepper_euler.do_step( lorenz, x1, t1, dt );
stepper_rk4.do_step( lorenz, x2, t2, dt );
file << " " << std::setw( 15 ) << std::setprecision( 10 ) << x1[0]
<< tab << std::setw( 15 ) << std::setprecision( 10 ) << x1[1]
<< tab << std::setw( 15 ) << std::setprecision( 10 ) << x1[2]
<< tab << "||"
<< " " << std::setw( 15 ) << std::setprecision( 10 ) << x2[0]
<< tab << std::setw( 15 ) << std::setprecision( 10 ) << x2[1]
<< tab << std::setw( 15 ) << std::setprecision( 10 ) << x2[2]
<< '\n';
}
file.close();
std::cout << "\nPress any key and enter to quit.\n";
std::cin.get();
return 0;
}
现在，您可以看到定义 Lorenz 函数并将其传递给您想要使用的 ODE 积分器类型是多么简单。
作为一个额外的好处，这里有一个很好的小函数，可以与这些结构一起使用来积分常量。
integrate_const.h
#ifndef INTEGRATE_CONST_H
#define INTEGRATE_CONST_H
// This function will iterate the state of the ODE,
// with the possibility to observe that state in each step
template<class Stepper, class DynamicSystem, class Observer>
size_t integrate_const( Stepper& stepper, DynamicSystem& system,
typename Stepper::container_type& state,
typename Stepper::time_type start_time,
typename Stepper::time_type end_time,
typename Stepper::time_type dt,
Observer& observer ) {
stepper.adjust_size( state );
size_t iteration = 0;
while( start_time < end_time ) {
observer( start_time, state, system );
stepper.do_step( system, state, start_time, dt );
start_time += dt;
++iteration;
}
observer( start_time, state, system );
return iteration;
}
#endif // !INTEGRATE_CONST_H
Summary- 现在可以使用任何具有开始、结束和调整大小功能的容器，而不是使用基本的数组或指针，从而使代码通用，易于实现和使用。上面的代码是工作代码，因为我没有看到任何编译或链接错误。我只对运行时错误进行了适度的测试。
Note:这项工作的灵感来自于我读过的一个项目code project。我确实对我找到的原始版本进行了一些修改，并且从boost图书馆。你可以找到它here https://www.codeproject.com/Articles/43607/Solving-ordinary-differential-equations-in-C。我也相信odeint库现在正式成为新版本的一部分boost图书馆。
以下是一些可能有用的有关使用常微分方程的其他链接：一些是关于对其进行编程的，另一些是关于其背后的数学符号的。
- 代码梦想 http://www.dreamincode.net/forums/topic/366744-runge-kutta-fourth-order-numerical-integrator/
- Github 上的 Odeint http://headmyshoulder.github.io/odeint-v2/
- 堆栈问答 https://stackoverflow.com/questions/42747227/c-numerical-integrators-to-solve-systems-of-odes
- You Tube https://www.youtube.com/user/HoustonMathPrep
