libode
Easy to compile, fast ODE integrators as C++ classes
Loading...
Searching...
No Matches
ode_vern_76.cc
Go to the documentation of this file.
1
2
3#include "ode_vern_76.h"
4
5namespace ode {
6
7OdeVern76::OdeVern76 (unsigned long neq) :
8 OdeEmbedded (neq, false, 6),
9 OdeRK (neq, 10),
10 OdeERK (neq) {
11
12 method_ = "Vern76";
13
14 /*
15 coefficients copied from:
16 http://people.math.sfu.ca/~jverner/
17 http://people.math.sfu.ca/~jverner/RKV76.IIa.Efficient.00001675585.081206.CoeffsOnlyFLOAT
18 */
19
20 c2 = 0.005;
21 c3 = 0.1088888888888888888888888888888888888889;
22 c4 = 0.1633333333333333333333333333333333333333;
23 c5 = 0.4555000000000000000000000000000000000000;
24 c6 = 0.6095094489978381317087004421486024949638;
25 c7 = 0.884;
26 c8 = 0.925;
27 c9 = 1.0;
28 c10 = 1.0;
29
30 a21 = c2;
31
32 a31 = -1.076790123456790123456790123456790123457;
33 a32 = 1.185679012345679012345679012345679012346;
34
35 a41 = 0.04083333333333333333333333333333333333333;
36 a43 = 0.1225;
37
38 a51 = 0.6389139236255726780508121615993336109954;
39 a53 = -2.455672638223656809662640566430653894211;
40 a54 = 2.272258714598084131611828404831320283215;
41
42 a61 = -2.661577375018757131119259297861818119279;
43 a63 = 10.80451388645613769565396655365532838482;
44 a64 = -8.353914657396199411968048547819291691541;
45 a65 = 0.8204875949566569791420417341743839209619;
46
47 a71 = 6.067741434696770992718360183877276714679;
48 a73 = -24.71127363591108579734203485290746001803;
49 a74 = 20.42751793078889394045773111748346612697;
50 a75 = -1.906157978816647150624096784352757010879;
51 a76 = 1.006172249242068014790040335899474187268;
52
53 a81 = 12.05467007625320299509109452892778311648;
54 a83 = -49.75478495046898932807257615331444758322;
55 a84 = 41.14288863860467663259698416710157354209;
56 a85 = -4.461760149974004185641911603484815375051;
57 a86 = 2.042334822239174959821717077708608543738;
58 a87 = -0.09834843665406107379530801693870224403537;
59
60 a91 = 10.13814652288180787641845141981689030769;
61 a93 = -42.64113603171750214622846006736635730625;
62 a94 = 35.76384003992257007135021178023160054034;
63 a95 = -4.348022840392907653340370296908245943710;
64 a96 = 2.009862268377035895441943593011827554771;
65 a97 = 0.3487490460338272405953822853053145879140;
66 a98 = -0.2714390051048312842371587140910297407572;
67
68 a101 = -45.03007203429867712435322405073769635151;
69 a103 = 187.3272437654588840752418206154201997384;
70 a104 = -154.0288236935018690596728621034510402582;
71 a105 = 18.56465306347536233859492332958439136765;
72 a106 = -7.141809679295078854925420496823551192821;
73 a107 = 1.308808578161378625114762706007696696508;
74
75 b1 = 0.04715561848627222170431765108838175679569;
76 b4 = 0.2575056429843415189596436101037687580986;
77 b5 = 0.2621665397741262047713863095764527711129;
78 b6 = 0.1521609265673855740323133199165117535523;
79 b7 = 0.4939969170032484246907175893227876844296;
80 b8 = -0.2943031171403250441557244744092703429139;
81 b9 = 0.08131747232495109999734599440136761892478;
82
83 d1 = 0.04460860660634117628731817597479197781432;
84 d4 = 0.2671640378571372680509102260943837899738;
85 d5 = 0.2201018300177293019979715776650753096323;
86 d6 = 0.2188431703143156830983120833512893824578;
87 d7 = 0.2289871705411202883378173889763552365362;
88 d10 = 0.02029518466335628222767054793810430358554;
89}
90
91void OdeVern76::step_ (double dt) {
92
93 unsigned long i;
94
95 //------------------------------------------------------------------
96 ode_fun_(sol_, k_[0]);
97
98 //------------------------------------------------------------------
99 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*a21*k_[0][i];
100 ode_fun_(soltemp_, k_[1]);
101
102 //------------------------------------------------------------------
103 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a31*k_[0][i]
104 + a32*k_[1][i]);
105 ode_fun_(soltemp_, k_[2]);
106
107 //------------------------------------------------------------------
108 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a41*k_[0][i]
109 + a43*k_[2][i]);
110 ode_fun_(soltemp_, k_[3]);
111
112 //------------------------------------------------------------------
113 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a51*k_[0][i]
114 + a53*k_[2][i]
115 + a54*k_[3][i]);
116 ode_fun_(soltemp_, k_[4]);
117
118 //------------------------------------------------------------------
119 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a61*k_[0][i]
120 + a63*k_[2][i]
121 + a64*k_[3][i]
122 + a65*k_[4][i]);
123 ode_fun_(soltemp_, k_[5]);
124
125 //------------------------------------------------------------------
126 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a71*k_[0][i]
127 + a73*k_[2][i]
128 + a74*k_[3][i]
129 + a75*k_[4][i]
130 + a76*k_[5][i]);
131 ode_fun_(soltemp_, k_[6]);
132
133 //------------------------------------------------------------------
134 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a81*k_[0][i]
135 + a83*k_[2][i]
136 + a84*k_[3][i]
137 + a85*k_[4][i]
138 + a86*k_[5][i]
139 + a87*k_[6][i]);
140 ode_fun_(soltemp_, k_[7]);
141
142 //------------------------------------------------------------------
143 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a91*k_[0][i]
144 + a93*k_[2][i]
145 + a94*k_[3][i]
146 + a95*k_[4][i]
147 + a96*k_[5][i]
148 + a97*k_[6][i]
149 + a98*k_[7][i]);
150 ode_fun_(soltemp_, k_[8]);
151
152 //------------------------------------------------------------------
153 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a101*k_[0][i]
154 + a103*k_[2][i]
155 + a104*k_[3][i]
156 + a105*k_[4][i]
157 + a106*k_[5][i]
158 + a107*k_[6][i]);
159 ode_fun_(soltemp_, k_[9]);
160
161 //------------------------------------------------------------------
162 for (i=0; i<neq_; i++) {
163 solemb_[i] = sol_[i] + dt*(d1*k_[0][i]
164 + d4*k_[3][i]
165 + d5*k_[4][i]
166 + d6*k_[5][i]
167 + d7*k_[6][i]
168 + d10*k_[9][i]);
169 sol_[i] = sol_[i] + dt*(b1*k_[0][i]
170 + b4*k_[3][i]
171 + b5*k_[4][i]
172 + b6*k_[5][i]
173 + b7*k_[6][i]
174 + b8*k_[7][i]
175 + b9*k_[8][i]);
176 }
177}
178
179} // 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::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::OdeVern76::OdeVern76
OdeVern76(unsigned long neq)
constructs
Definition ode_vern_76.cc:7
ode_vern_76.h