libpspm
Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
RKCK45 Class Reference

#include <rkck45.h>

Collaboration diagram for RKCK45:
[legend]

Public Member Functions

 RKCK45 (double t_start_, double accuracy, double h1)
 
void resize (int new_size)
 
template<class functor >
void Step (double &x, container &y, functor &derivs, double hmax=1e20, double hmin=1e-6)
 
template<class functor , class AfterStep >
void Step_to (double t_stop, double &x, container &y, functor &derivs, AfterStep &after_step)
 
double h ()
 
double size ()
 
int get_fn_evals ()
 

Private Member Functions

template<class functor >
void RKTry (container &y, container &dydx, double &x, double h, container &yout, container &yerr, functor &derivs)
 
template<class functor >
void RKStep (container &y, container &dydx, double &x, double htry, double &hdid, double &hnext, functor &derivs)
 

Private Attributes

double a_y = 1
 
double a_dydt = 0
 
container yscal
 
double ht
 
double eps_rel
 
double eps_abs
 
double xt
 
double t_stop
 
int nok =0
 
int nbad =0
 
container dydx
 
int nfe = 0
 
container k1
 
container k2
 
container k3
 
container k4
 
container k5
 
container yt
 
int sys_size = 0
 

Static Private Attributes

static constexpr double SAFETY = 0.9
 
static constexpr double PGROW = -0.2
 
static constexpr double PSHRNK = -0.25
 
static constexpr double ERRCON = 1.89e-4
 

Detailed Description

Definition at line 47 of file rkck45.h.

Constructor & Destructor Documentation

◆ RKCK45()

RKCK45::RKCK45 ( double  t_start_,
double  accuracy,
double  h1 
)

Definition at line 5 of file rkck45.cpp.

5  :
6  ht(h1), eps_rel(accuracy), eps_abs(accuracy), xt(t_start_){
7 
8 }
RKCK45::eps_rel
double eps_rel
Definition: rkck45.h:57
RKCK45::ht
double ht
Definition: rkck45.h:56
RKCK45::eps_abs
double eps_abs
Definition: rkck45.h:57
RKCK45::xt
double xt
Definition: rkck45.h:58

Member Function Documentation

◆ get_fn_evals()

int RKCK45::get_fn_evals ( )
inline

Definition at line 96 of file rkck45.h.

96 {return nfe;}
RKCK45::nfe
int nfe
Definition: rkck45.h:61
Here is the call graph for this function:

◆ h()

double RKCK45::h ( )
inline

Definition at line 94 of file rkck45.h.

94 {return ht;}
RKCK45::ht
double ht
Definition: rkck45.h:56
Here is the caller graph for this function:

◆ resize()

void RKCK45::resize ( int  new_size)

Definition at line 10 of file rkck45.cpp.

10  {
11  sys_size = new_size;
12 
13  yscal.resize(new_size);
14  dydx.resize(new_size);
15  k1.resize(new_size);
16  k2.resize(new_size);
17  k3.resize(new_size);
18  k4.resize(new_size);
19  k5.resize(new_size);
20  yt.resize(new_size);
21 }
RKCK45::yt
container yt
Definition: rkck45.h:63
RKCK45::k3
container k3
Definition: rkck45.h:63
RKCK45::yscal
container yscal
Definition: rkck45.h:55
RKCK45::k4
container k4
Definition: rkck45.h:63
RKCK45::k1
container k1
Definition: rkck45.h:63
RKCK45::k2
container k2
Definition: rkck45.h:63
RKCK45::k5
container k5
Definition: rkck45.h:63
RKCK45::dydx
container dydx
Definition: rkck45.h:60
RKCK45::sys_size
int sys_size
Definition: rkck45.h:65
Here is the caller graph for this function:

◆ RKStep()

template<class functor >
void RKCK45::RKStep ( container &  y,
container &  dydx,
double &  x,
double  htry,
double &  hdid,
double &  hnext,
functor &  derivs 
)
private

Definition at line 11 of file rkck45.tpp.

12  {
13 
14  container yerr(y.size()), ytemp(y.size()); // TODO: Should these also be made class members to prevent reallocation?
15  double errmax;
16  double h=htry; // Set stepsize to the initial trial value.
17  for (;;) { // infinite loop
18  RKTry(y, dydx, x, h, ytemp, yerr, derivs);// Take a step.
19 
20  errmax=0.0; // Evaluate accuracy.
21  for (int i=0; i<y.size(); i++) errmax=std::max(errmax, fabs(yerr[i]/yscal[i]));
22  //errmax /= eps; // Scale relative to required tolerance.
23  //std:: cout << "x = " << x << ", errmax = " << errmax << "\n";
24  if (errmax <= 1.1) break; // Step succeeded. Compute size of next step.
25  double htemp=h*fmax(SAFETY*pow(errmax,PSHRNK), 0.2); // Truncation error too large, reduce stepsize, max by a factor of 0.2.
26  h=htemp; //(h >= 0.0 ? std::max(htemp,0.1*h) : std::min(htemp,0.1*h)); // No more than a factor of 10.
27  if ((x+h) == x) std::cerr<<"stepsize underflow in rkqs"<<std::endl;
28  }
29  if (errmax < 0.5) hnext=h*fmin(SAFETY*pow(errmax,PGROW), 5); // Step is too small, increase it next time, no more than 5 times
30  else hnext=h;
31  //cout << "hnext = " << hnext << "\n";
32  hdid=h;
33  x += h;
34  for (int i=0; i<y.size();i++) y[i] = ytemp[i];
35 }
RKCK45::PGROW
static constexpr double PGROW
Definition: rkck45.h:50
RKCK45::yscal
container yscal
Definition: rkck45.h:55
RKCK45::h
double h()
Definition: rkck45.h:94
RKCK45::PSHRNK
static constexpr double PSHRNK
Definition: rkck45.h:51
RKCK45::SAFETY
static constexpr double SAFETY
Definition: rkck45.h:49
RKCK45::RKTry
void RKTry(container &y, container &dydx, double &x, double h, container &yout, container &yerr, functor &derivs)
Definition: rkck45.tpp:43
RKCK45::dydx
container dydx
Definition: rkck45.h:60
container
std::vector< double > container
Definition: rkck45.h:32
Here is the call graph for this function:
Here is the caller graph for this function:

◆ RKTry()

template<class functor >
void RKCK45::RKTry ( container &  y,
container &  dydx,
double &  x,
double  h,
container &  yout,
container &  yerr,
functor &  derivs 
)
private

Definition at line 43 of file rkck45.tpp.

43  {
44  //std::cout << "try:\n";
45  // a0 a1 a2 a3 a4 a5
46  static constexpr double as[] = {0, 0.2, 0.3, 0.6, 1.0, 0.875};
47  // c0 c1 c2 c3 c4 c5
48  static constexpr double cs[] = {37.0/378.0, 0.0, 250.0/621.0, 125.0/594.0, 0.0, 512.0/1771.0};
49  // dc0 dc1 dc2 dc3 dc4 dc5
50  static constexpr double dc[] = {cs[0]-2825.0/27648.0, 0.0, cs[2]-18575.0/48384.0, cs[3]-13525.0/55296.0, -277.00/14336.0, cs[5]-0.25};
51  static constexpr double bs[6][6] =
52  {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
53  {0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
54  {3.0/40.0, 9.0/40.0, 0.0, 0.0, 0.0, 0.0},
55  {0.3, -0.9, 1.2, 0.0, 0.0, 0.0},
56  {-11.0/54.0, 2.5, -70.0/27.0, 35.0/27.0, 0.0, 0.0},
57  {1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0, 0.0}};
58 
59  int N=y.size();
60 
61  container& k0 = dydx; // just different name for the same variable
62 
63  // First step.
64  for (int i=0; i<N; i++) yt[i] = y[i] + h*bs[1][0]*k0[i];
65  derivs(x+as[1]*h, yt.begin(), k1.begin(), nullptr);
66  ++nfe;
67 
68  // Second step.
69  for (int i=0; i<N; i++) yt[i] = y[i] + h*(bs[2][0]*k0[i]+bs[2][1]*k1[i]);
70  derivs(x+as[2]*h,yt.begin(),k2.begin(), nullptr);
71  ++nfe;
72 
73  // Third step.
74  for (int i=0; i<N; i++) yt[i] = y[i] + h*(bs[3][0]*k0[i]+bs[3][1]*k1[i]+bs[3][2]*k2[i]);
75  derivs(x+as[3]*h,yt.begin(),k3.begin(), nullptr);
76  ++nfe;
77 
78  // Fourth step.
79  for (int i=0; i<N; i++) yt[i] = y[i] + h*(bs[4][0]*k0[i]+bs[4][1]*k1[i]+bs[4][2]*k2[i]+bs[4][3]*k3[i]);
80  derivs(x+as[4]*h,yt.begin(),k4.begin(), nullptr);
81  ++nfe;
82 
83  // Fifth step.
84  for (int i=0;i<N;i++) yt[i] = y[i] + h*(bs[5][0]*k0[i]+bs[5][1]*k1[i]+bs[5][2]*k2[i]+bs[5][3]*k3[i]+bs[5][4]*k4[i]);
85  derivs(x+as[5]*h,yt.begin(),k5.begin(), nullptr);
86  ++nfe;
87 
88  // Sixth step.
89  // Accumulate increments with proper weights.
90  for (int i=0;i<N;i++) yout[i] = y[i] + h*(cs[0]*k0[i]+cs[2]*k2[i]+cs[3]*k3[i]+cs[5]*k5[i]);
91 
92  // Estimate error as difference between fourth and fifth order methods.
93  for (int i=0; i<N; i++) yerr[i] = h*(dc[0]*dydx[i]+dc[2]*k2[i]+dc[3]*k3[i]+dc[4]*k4[i]+dc[5]*k5[i]);
94 }
RKCK45::yt
container yt
Definition: rkck45.h:63
RKCK45::k3
container k3
Definition: rkck45.h:63
RKCK45::h
double h()
Definition: rkck45.h:94
RKCK45::nfe
int nfe
Definition: rkck45.h:61
RKCK45::k4
container k4
Definition: rkck45.h:63
RKCK45::k1
container k1
Definition: rkck45.h:63
RKCK45::k2
container k2
Definition: rkck45.h:63
RKCK45::k5
container k5
Definition: rkck45.h:63
RKCK45::dydx
container dydx
Definition: rkck45.h:60
container
std::vector< double > container
Definition: rkck45.h:32
Here is the caller graph for this function:

◆ size()

double RKCK45::size ( )
inline

Definition at line 95 of file rkck45.h.

95 {return sys_size;}
RKCK45::sys_size
int sys_size
Definition: rkck45.h:65

◆ Step()

template<class functor >
void RKCK45::Step ( double &  x,
container &  y,
functor &  derivs,
double  hmax = 1e20,
double  hmin = 1e-6 
)

Definition at line 103 of file rkck45.tpp.

103  {
104 
105  // resize the container if necessary
106  if (y.size() != sys_size) resize(y.size());
107 
108  // Calculates derivatives for the new step
109  //cout << "\n>> init derivs ~~~\n";
110  derivs(xt,y.begin(),dydx.begin(), nullptr);
111  ++nfe;
112  // good way of determining desired accuracy
113  for (int i=0; i<y.size(); i++) yscal[i] = eps_abs + eps_rel*(a_y * fabs(y[i]) + a_dydt * fabs(dydx[i]*ht)); // fabs(y[i]) + fabs(dydx[i]*ht) + 1e-3;
114  double hnext, hdid;
115  // Cals the Runge-Kutta rountine
116  // takes: y -- dependent variable, dydx -- derivative at the beginning
117  // ht -- trial setpsize, hdid --
118  ht = std::min(ht, hmax);
119  ht = std::max(ht, hmin);
120  RKStep(y, dydx, xt, ht, hdid, hnext, derivs); // Make one RK step
121  if (hdid == ht) ++nok; else ++nbad; // good or bad step
122  ht = hnext;
123  x = xt;
124  //return x<t_stop;
125 }
RKCK45::nok
int nok
Definition: rkck45.h:59
RKCK45::a_dydt
double a_dydt
Definition: rkck45.h:54
RKCK45::eps_rel
double eps_rel
Definition: rkck45.h:57
RKCK45::nbad
int nbad
Definition: rkck45.h:59
RKCK45::yscal
container yscal
Definition: rkck45.h:55
RKCK45::RKStep
void RKStep(container &y, container &dydx, double &x, double htry, double &hdid, double &hnext, functor &derivs)
Definition: rkck45.tpp:11
RKCK45::nfe
int nfe
Definition: rkck45.h:61
RKCK45::ht
double ht
Definition: rkck45.h:56
RKCK45::resize
void resize(int new_size)
Definition: rkck45.cpp:10
RKCK45::eps_abs
double eps_abs
Definition: rkck45.h:57
RKCK45::xt
double xt
Definition: rkck45.h:58
RKCK45::a_y
double a_y
Definition: rkck45.h:53
RKCK45::dydx
container dydx
Definition: rkck45.h:60
RKCK45::sys_size
int sys_size
Definition: rkck45.h:65
Here is the call graph for this function:
Here is the caller graph for this function:

◆ Step_to()

template<class functor , class AfterStep >
void RKCK45::Step_to ( double  t_stop,
double &  x,
container &  y,
functor &  derivs,
AfterStep &  after_step 
)

Definition at line 149 of file rkck45.tpp.

149  {
150  while (x < t_stop){
151  Step(x,y,derivs, t_stop-x);
152  after_step(x, y.begin());
153  }
154 }
RKCK45::Step
void Step(double &x, container &y, functor &derivs, double hmax=1e20, double hmin=1e-6)
Definition: rkck45.tpp:103
RKCK45::t_stop
double t_stop
Definition: rkck45.h:58
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ a_dydt

double RKCK45::a_dydt = 0
private

Definition at line 54 of file rkck45.h.

◆ a_y

double RKCK45::a_y = 1
private

Definition at line 53 of file rkck45.h.

◆ dydx

container RKCK45::dydx
private

Definition at line 60 of file rkck45.h.

◆ eps_abs

double RKCK45::eps_abs
private

Definition at line 57 of file rkck45.h.

◆ eps_rel

double RKCK45::eps_rel
private

Definition at line 57 of file rkck45.h.

◆ ERRCON

constexpr double RKCK45::ERRCON = 1.89e-4
staticprivate

Definition at line 52 of file rkck45.h.

◆ ht

double RKCK45::ht
private

Definition at line 56 of file rkck45.h.

◆ k1

container RKCK45::k1
private

Definition at line 63 of file rkck45.h.

◆ k2

container RKCK45::k2
private

Definition at line 63 of file rkck45.h.

◆ k3

container RKCK45::k3
private

Definition at line 63 of file rkck45.h.

◆ k4

container RKCK45::k4
private

Definition at line 63 of file rkck45.h.

◆ k5

container RKCK45::k5
private

Definition at line 63 of file rkck45.h.

◆ nbad

int RKCK45::nbad =0
private

Definition at line 59 of file rkck45.h.

◆ nfe

int RKCK45::nfe = 0
private

Definition at line 61 of file rkck45.h.

◆ nok

int RKCK45::nok =0
private

Definition at line 59 of file rkck45.h.

◆ PGROW

constexpr double RKCK45::PGROW = -0.2
staticprivate

Definition at line 50 of file rkck45.h.

◆ PSHRNK

constexpr double RKCK45::PSHRNK = -0.25
staticprivate

Definition at line 51 of file rkck45.h.

◆ SAFETY

constexpr double RKCK45::SAFETY = 0.9
staticprivate

Definition at line 49 of file rkck45.h.

◆ sys_size

int RKCK45::sys_size = 0
private

Definition at line 65 of file rkck45.h.

◆ t_stop

double RKCK45::t_stop
private

Definition at line 58 of file rkck45.h.

◆ xt

double RKCK45::xt
private

Definition at line 58 of file rkck45.h.

◆ yscal

container RKCK45::yscal
private

Definition at line 55 of file rkck45.h.

◆ yt

container RKCK45::yt
private

Definition at line 63 of file rkck45.h.

The documentation for this class was generated from the following files: