在C ++中使用RK-4求解Lorenz方程 [英] Solving Lorenz equation using RK-4 in C++
本文介绍了在C ++中使用RK-4求解Lorenz方程的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
我编写了使用C ++中的RK-4方法求解Lorenz方程的代码.我必须绘制吸引子图，并且在使用RK-4方法求解3个一阶耦合微分方程时遇到一些困难.这是我的代码:
I wrote the code to solve Lorenz equations using RK-4 method in C++. I have to plot the attractor plot and I have some difficulty in solving the 3 first order coupled differential equation using RK-4 method. Here is my code:
/* 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);
}
编译时，我收到一条错误消息，说该程序停止工作.有人可以帮我吗?
When I compile it, I get an error saying that program stopped working. Can someone please help me out?
推荐答案
我问过OP他们是否能够使用模板，他们回答"I can use templates, but I wanted it to be simple."模板结构；但是一旦完成；模板的使用确实简化了很多事情.最后，它确实使代码更通用，而不是特定于单个任务.
I had asked the OP if they were able to use templates or not, they replied with "I can use templates, but I wanted it to be simple." There is a little more work behind the scenes to first set up the template structures; but once that is done; the use of templates does simplify things quite a bit. In the end it does generally make the code more generic as opposed to being specific to a single task.
这是一个未实现的类，用于显示ODE集成器的通用模式.
Here is a non implemented class to show the generalized pattern of an ODE integrator.
#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进行建模之前，我们需要定义一些东西.
The above will not do anything but only serves to show how an ODE is to behave or act. Before we can model any kind of ODE after the above pattern we do need to define a few things.
首先，我们需要一个container_traits类.
First we need a container_traits class.
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
上面的类很简单，很简单.
The above class is fairly simple and pretty straight forward.
然后，我们需要一组迭代器类型的函数来迭代那些容器.迭代器的功能要稍微复杂一些，但是对于它们的预期目的，它还是很简单的.我将它们包装在名称空间中，这样它们就不会与可能存在于其中的任何其他迭代器类或函数发生冲突.
Then we need a set of iterator type functions to iterate those containers. The iterator functions are a little more involved but again it is fairly straight forward as to what their intended purposes are. I have them wrapped in a namespace so they won't conflict with any other iterator classes or functions that may be out there.
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或其他任何形式.
Now that we have our needed structure's and functions we can use the pattern above to implement different types of ODE integrators. I will show two of the most common: euler and rk4. One can easily implement midpoint, rk5 or any other as long as they follow the pattern.
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
现在，我们拥有了所需的一切；剩下要做的就是使用它们.
Now that we have all that we need; all that is left to do is to use them.
OP提到要解决Lorenz问题，因此我将用它来演示如何使用ODE:
The OP mentioned about solving a Lorenz problem so I will use that to demonstrate how to use the ODEs:
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集成器的类型是多么简单.
You can now see how simple it is to define a Lorenz function and simply pass that to the type of ODE integrator you want to use.
作为一个额外的好处，这是一个很好的小函数，可以与这些结构一起使用以集成常量.
As an added bonus here is a nice little function that can be used with these structures for integrating constants.
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 - Instead of using basic arrays or pointers, one can now use any container that has a begin, end & resize function making the code generic, easy to implement and use. The above code is working code as I have not seen any compile nor linking errors. I have only moderately tested for runtime errors. If anyone happens to find any bugs; please do not down vote because of a bug, but feel free to leave a comment describing the bug you found and I will gladly fix the appropriate bug(s) to make this code base more accurate.
注意:这项工作的灵感来自于我从code project中读取的一个项目.我确实对找到的原始版本进行了一些修改，并从boost库中删除了所有依赖项.您可以在此处找到它.我也相信odeint库现在正式是boost库的较新版本的一部分.
Note: This work was inspired by a project that I read from code project. I did make a few modifications from the original version that I've found, and I removed any dependencies from the boost library. You can find it here. Also I believe that the odeint library is now officially a part of the newer versions of the boost library.
这里还有一些其他有关使用ODE的链接，这些链接可能会有用:一些与编程有关，而另一些与背后的数学符号有关.
Here are some other links about working with ODEs that may be useful: Some are about programming it others are about the math notation behind it.
- Dream In Code
- You Tube
- Odeint on Github
- Stack Q/A
- You Tube
这篇关于在C ++中使用RK-4求解Lorenz方程的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
