我编写了代码，用C ++中的RK-4方法解决Lorenz方程。我必须绘制吸引子图，并且在使用RK-4方法求解3阶一阶耦合微分方程时遇到一些困难。这是我的代码：
/* Solving 3 coupled first order differential equations using RK4 for the
Lorentz system
dxdt = sig*(y - x) {f}
dydt = x*(rho - z) - y {g}
dzdt = x*y - bet*z {k} */
#include<iostream>
#include<cmath>
#include<fstream>
#include <iomanip>
using namespace std;
double sig = 10.0;
double rho = 28.0;
double bet = 8.0/3.0;
double func1 (double y[], int N) // function to define the n first order
differential equations
{
double dxdt;
dxdt = sig*(y[1] - y[0]);
return dxdt;
}
double func2 (double y[], int N) // function to define the n first order
differential equations
{
double dydt;
dydt = y[0]*(rho - y[2]) - y[1];
return dydt;
}
double func3 (double y[], int N) // function to define the n first order
differential equations
{
double dzdt;
dzdt = y[0]*y[1] - bet*y[2];
return dzdt;
}
void rk4(int n, double x0,double xf,double Y0[],int N) // Function to
implement RK4 algorithm
{
double K1;
double K2;
double K3;
double K4;
double L1;
double L2;
double L3;
double L4;
double M1;
double M2;
double M3;
double M4;
double x[n+1];
double h = (xf-x0)/n;
long double y[n+1][N];
double dydx[N];
double ytemp[N];
for(int i=0;i<=n;i++) // loop to store values of time
{
x[i] = x0 + i*h;
}
for(int j =0;j<N;j++) // loop to store initial conditions of diff eq
{
y[0][j] = Y0[j];
}
for (int k =0;k<n;k++)
{
for(int l=0;l<N;l++)
{
ytemp[l] = y[k][l];
}
K1 = func1(ytemp,N); // f
L1 = func2(ytemp,N); // g
M1 = func3(ytemp,N); // k
ytemp[0] = y[k][0] + 0.5*h*K1;
ytemp[1] = y[k][1] + 0.5*h*L1;
ytemp[2] = y[k][2] + 0.5*h*M1;
K2 = func1(ytemp,N);
L2 = func2(ytemp,N);
M2 = func3(ytemp,N);
ytemp[0] = y[k][0] + 0.5*h*K2;
ytemp[1] = y[k][1] + 0.5*h*L2;
ytemp[2] = y[k][2] + 0.5*h*M2;
K3 = func1(ytemp,N);
L3 = func2(ytemp,N);
M3 = func3(ytemp,N);
ytemp[0] = y[k][0] + h*K3;
ytemp[1] = y[k][1] + h*L3;
ytemp[2] = y[k][2] + h*M3;
K4 = func1(ytemp,N);
L4 = func2(ytemp,N);
M4 = func3(ytemp,N);
dydx[0] = (1.0/6.0)*(K1+ 2.0*K2 + 2.0*K3+ K4);
dydx[1] = (1.0/6.0)*(L1+ 2.0*L2 + 2.0*L3+ L4);
dydx[2] = (1.0/6.0)*(M1+ 2.0*M2 + 2.0*M3+ M4);
for(int l=0;l<N;l++)
{
y[k+1][l] = y[k][l] + h*dydx[l];
}
}
ofstream fout;
fout.open("prog12bdata.txt",ios::out);
for (int m =0;m<=n;m++)
{
fout<<setw(15) << setprecision(10) <<x[m];
for(int o =0;o<N;o++)
{
fout<<" "<<setw(15) << setprecision(10) <<y[m][o];
}
fout<<endl;
}
fout.close();
}
int main()
{
int N;// Order of ODE to solve
cout<<"Enter the order of ODE to solve: "<<endl;
cin>>N;
double x0=0;
double xf=50;
int n = 50000; // number of steps
double Y0[N];
for(int i=0;i<N;i++)
{
cout<<"Enter the initial conditions: "<<endl;
cin>>Y0[i];
}
rk4(n,x0,xf,Y0,N);
}
当我编译它时，我得到一个错误，说程序停止工作。有人可以帮帮我吗？
我问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。只要符合模式，就可以轻松实现中点，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提到了解决Lorenz问题所以我将用它来演示如何使用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
总结 - 现在可以使用任何具有begin，end和resize函数的容器，而不是使用基本数组或指针，使代码具有通用性，易于实现和使用。上面的代码是工作代码，因为我没有看到任何编译或链接错误。我只对运行时错误进行了适度测试。如果有人碰巧发现任何错误;请不要因为一个错误而投票，但随时留下评论描述你发现的错误，我很乐意修复相应的错误，使这个代码库更准确。
注意：这项工作的灵感来自我从code project读到的一个项目。我确实从我发现的原始版本做了一些修改，并从boost库中删除了任何依赖项。你可以找到它here。另外我相信odeint库现在正式成为boost库的新版本的一部分。
以下是一些关于使用可能有用的ODE的其他链接：一些是关于编程的，另一些是关于它背后的数学符号。