libode
Easy to compile, fast ODE integrators as C++ classes
Loading...
Searching...
No Matches
ode_dopri_87.cc
Go to the documentation of this file.
1
2
3#include "ode_dopri_87.h"
4
5namespace ode {
6
7OdeDoPri87::OdeDoPri87 (unsigned long neq) :
8 OdeEmbedded (neq, false, 7),
9 OdeRK (neq, 13),
10 OdeERK (neq) {
11
12 method_ = "DoPri87";
13 //tableau of coefficents
14 c2 = 1.0/18; a21 = c2;
15 c3 = 1.0/12; a31 = 1.0/48; a32 = 1.0/16;
16 c4 = 1.0/8; a41 = 1.0/32; a43 = 3.0/32;
17 c5 = 5.0/16; a51 = 5.0/16; a53 = -75.0/64; a54 = 75.0/64;
18 c6 = 3.0/8; a61 = 3.0/80; a64 = 3.0/16; a65 = 3.0/20;
19 c7 = 59.0/400; a71 = 29443841.0/614563906; a74 = 77736538.0/692538347; a75 = -28693883.0/1125000000; a76 = 23124283.0/1800000000;
20 c8 = 93.0/200; a81 = 16016141.0/946692911; a84 = 61564180.0/158732637; a85 = 22789713.0/633445777; a86 = 545815736.0/2771057229; a87 = -180193667.0/1043307555;
21 c9 = 5490023248.0/9719169821; a91 = 39632708.0/573591083; a94 = -433636366.0/683701615; a95 = -421739975.0/2616292301; a96 = 100302831.0/723423059; a97 = 790204164.0/839813087; a98 = 800635310.0/3783071287;
22 c10 = 13.0/20; a101 = 246121993.0/1340847787; a104 = -37695042795.0/15268766246; a105 = -309121744.0/1061227803; a106 = -12992083.0/490766935; a107 = 6005943493.0/2108947869; a108 = 393006217.0/1396673457; a109 = 123872331.0/1001029789;
23 c11 = 1201146811.0/1299019798; a111 = -1028468189.0/846180014; a114 = 8478235783.0/508512852; a115 = 1311729495.0/1432422823; a116 = -10304129995.0/1701304382; a117 = -48777925059.0/3047939560; a118 = 15336726248.0/1032824649; a119 = -45442868181.0/3398467696; a1110 = 3065993473.0/597172653;
24 c12 = 1.0; a121 = 185892177.0/718116043; a124 = -3185094517.0/667107341; a125 = -477755414.0/1098053517; a126 = -703635378.0/230739211; a127 = 5731566787.0/1027545527; a128 = 5232866602.0/850066563; a129 = -4093664535.0/808688257; a1210 = 3962137247.0/1805957418; a1211 = 65686358.0/487910083;
25 c13 = 1.0; a131 = 403863854.0/491063109; a134 = -5068492393.0/434740067; a135 = -411421997.0/543043805; a136 = 652783627.0/914296604; a137 = 11173962825.0/925320556; a138 = -13158990841.0/6184727034; a139 = 3936647629.0/1978049680; a1310 = -160528059.0/685178525; a1311 = 248638103.0/1413531060;
26 b1 = 14005451.0/335480064; b6 = -59238493.0/1068277825; b7 = 181606767.0/758867731; b8 = 561292985.0/797845732; b9 = -1041891430.0/1371343529; b10 = 760417239.0/1151165299; b11 = 118820643.0/751138087; b12 = -528747749.0/2220607170; b13 = 1.0/4;
27 d1 = 13451932.0/455176632; d6 = -808719846.0/976000145; d7 = 1757004468.0/5645159321; d8 = 656045339.0/265891186; d9 = -3867574721.0/1518517206; d10 = 465885868.0/322736535; d11 = 53011238.0/667516719; d12 = 2.0/45;
28}
29
30void OdeDoPri87::step_ (double dt) {
31
32 unsigned long i;
33
34 //------------------------------------------------------------------
35 ode_fun_(sol_, k_[0]);
36
37 //------------------------------------------------------------------
38 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*a21*k_[0][i];
39 ode_fun_(soltemp_, k_[1]);
40
41 //------------------------------------------------------------------
42 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a31*k_[0][i]
43 + a32*k_[1][i]);
44 ode_fun_(soltemp_, k_[2]);
45
46 //------------------------------------------------------------------
47 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a41*k_[0][i]
48 + a43*k_[2][i]);
49 ode_fun_(soltemp_, k_[3]);
50
51 //------------------------------------------------------------------
52 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a51*k_[0][i]
53 + a53*k_[2][i]
54 + a54*k_[3][i]);
55 ode_fun_(soltemp_, k_[4]);
56
57 //------------------------------------------------------------------
58 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a61*k_[0][i]
59 + a64*k_[3][i]
60 + a65*k_[4][i]);
61 ode_fun_(soltemp_, k_[5]);
62
63 //------------------------------------------------------------------
64 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a71*k_[0][i]
65 + a74*k_[3][i]
66 + a75*k_[4][i]
67 + a76*k_[5][i]);
68 ode_fun_(soltemp_, k_[6]);
69
70 //------------------------------------------------------------------
71 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a81*k_[0][i]
72 + a84*k_[3][i]
73 + a85*k_[4][i]
74 + a86*k_[5][i]
75 + a87*k_[6][i]);
76 ode_fun_(soltemp_, k_[7]);
77
78 //------------------------------------------------------------------
79 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a91*k_[0][i]
80 + a94*k_[3][i]
81 + a95*k_[4][i]
82 + a96*k_[5][i]
83 + a97*k_[6][i]
84 + a98*k_[7][i]);
85 ode_fun_(soltemp_, k_[8]);
86
87 //------------------------------------------------------------------
88 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a101*k_[0][i]
89 + a104*k_[3][i]
90 + a105*k_[4][i]
91 + a106*k_[5][i]
92 + a107*k_[6][i]
93 + a108*k_[7][i]
94 + a109*k_[8][i]);
95 ode_fun_(soltemp_, k_[9]);
96
97 //------------------------------------------------------------------
98 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a111*k_[0][i]
99 + a114*k_[3][i]
100 + a115*k_[4][i]
101 + a116*k_[5][i]
102 + a117*k_[6][i]
103 + a118*k_[7][i]
104 + a119*k_[8][i]
105 + a1110*k_[9][i]);
106 ode_fun_(soltemp_, k_[10]);
107
108 //------------------------------------------------------------------
109 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a121*k_[0][i]
110 + a124*k_[3][i]
111 + a125*k_[4][i]
112 + a126*k_[5][i]
113 + a127*k_[6][i]
114 + a128*k_[7][i]
115 + a129*k_[8][i]
116 + a1210*k_[9][i]
117 + a1211*k_[10][i]);
118 ode_fun_(soltemp_, k_[11]);
119
120 //------------------------------------------------------------------
121 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a131*k_[0][i]
122 + a134*k_[3][i]
123 + a135*k_[4][i]
124 + a136*k_[5][i]
125 + a137*k_[6][i]
126 + a138*k_[7][i]
127 + a139*k_[8][i]
128 + a1310*k_[9][i]
129 + a1311*k_[10][i]);
130 ode_fun_(soltemp_, k_[12]);
131
132 //------------------------------------------------------------------
133 for (i=0; i<neq_; i++) {
134 solemb_[i] = sol_[i] + dt*(d1*k_[0][i]
135 + d6*k_[5][i]
136 + d7*k_[6][i]
137 + d8*k_[7][i]
138 + d9*k_[8][i]
139 + d10*k_[9][i]
140 + d11*k_[10][i]
141 + d12*k_[11][i]);
142 sol_[i] = sol_[i] + dt*(b1*k_[0][i]
143 + b6*k_[5][i]
144 + b7*k_[6][i]
145 + b8*k_[7][i]
146 + b9*k_[8][i]
147 + b10*k_[9][i]
148 + b11*k_[10][i]
149 + b12*k_[11][i]
150 + b13*k_[12][i]);
151 }
152}
153
154} // namespace ode
ode::OdeBase::neq_
unsigned long neq_
number of equations in the system of ODEs
Definition ode_base.h:327
ode::OdeBase::method_
std::string method_
the "name" of the solver/method, as in "Euler" or "RK4"
Definition ode_base.h:319
ode::OdeBase::sol_
double * sol_
array for the solution, changing over time
Definition ode_base.h:333
ode::OdeBase::ode_fun_
void ode_fun_(double *solin, double *fout)
wrapper, calls ode_fun() and increases the neval counter by one
Definition ode_base.cc:87
ode::OdeDoPri87::OdeDoPri87
OdeDoPri87(unsigned long neq)
constructs
Definition ode_dopri_87.cc:7
ode::OdeERK
Base class providing space for temporary solutions moving through RK stages.
Definition ode_erk.h:9
ode::OdeERK::soltemp_
double * soltemp_
temporary solution vector
Definition ode_erk.h:22
ode::OdeEmbedded
Base clase implementing methods for embedded Runge-Kutta error estimation.
Definition ode_embedded.h:15
ode::OdeEmbedded::solemb_
double * solemb_
embedded solution array
Definition ode_embedded.h:63
ode::OdeRK
Provides space for stage slope values, an array of arrays for k values.
Definition ode_rk.h:9
ode::OdeRK::k_
double ** k_
stage evaluations
Definition ode_rk.h:23
ode_dopri_87.h