1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
/*
* rt_explicit_rk.hpp
*
* Copyright 2011 Mario Mulansky
* Copyright 2012 Karsten Ahnert
*
* Distributed under the Boost Software License, Version 1.0.
* (See accompanying file LICENSE_1_0.txt or
* copy at http://www.boost.org/LICENSE_1_0.txt)
*/
#ifndef RT_EXPLICIT_RK_HPP_
#define RT_EXPLICIT_RK_HPP_
#include <vector>
#include "rt_algebra.hpp"
using namespace std;
template< class StateType >
class rt_explicit_rk
{
public:
typedef StateType state_type;
typedef double* const * coeff_a_type;
typedef vector< double > coeff_b_type;
typedef vector< double > coeff_c_type;
rt_explicit_rk( size_t stage_count ) : m_s( stage_count )
{
m_F = new state_type[ m_s ];
}
rt_explicit_rk( const size_t stage_count ,
const coeff_a_type a ,
const coeff_b_type &b , const coeff_c_type &c )
: m_s( stage_count ) , m_a( a ) , m_b( b ) , m_c( c )
{
m_F = new state_type[ m_s ];
}
~rt_explicit_rk()
{
delete[] m_F;
}
/* void set_params( coeff_a_type &a , coeff_b_type &b , coeff_c_type &c )
{
m_a = a;
m_b = b;
m_c = c;
}*/
template< class System >
void do_step( System sys , state_type &x , const double t , const double dt )
{
// first stage separately
sys( x , m_F[0] , t + m_c[0]*t );
if( m_s == 1 )
rt_algebra::foreach( x , x , &m_b[0] , m_F , dt , 1 );
else
rt_algebra::foreach( m_x_tmp , x , m_a[0] , m_F , dt , 1 );
for( size_t stage = 2 ; stage <= m_s ; ++stage )
{
sys( m_x_tmp , m_F[stage-1] , t + m_c[stage-1]*dt );
if( stage == m_s )
rt_algebra::foreach( x , x , &m_b[0] , m_F , dt , stage-1 );
else
rt_algebra::foreach( m_x_tmp , x , m_a[stage-1] , m_F , dt , stage-1 );
}
}
private:
const size_t m_s;
const coeff_a_type m_a;
const coeff_b_type m_b;
const coeff_c_type m_c;
state_type m_x_tmp;
state_type *m_F;
};
#endif /* RT_EXPLICIT_RK_HPP_ */
|
