libode
Easy to compile, fast ODE integrators as C++ classes
Loading...
Searching...
No Matches
ode_vern_98.cc
Go to the documentation of this file.
1
2
3#include "ode_vern_98.h"
4
5namespace ode {
6
7OdeVern98::OdeVern98 (unsigned long neq) :
8 OdeEmbedded (neq, false, 8),
9 OdeRK (neq, 16),
10 OdeERK (neq) {
11
12 method_ = "Vern98";
13
14 /*
15 coefficients copied from:
16 http://people.math.sfu.ca/~jverner/
17 http://people.math.sfu.ca/~jverner/RKV98.IIa.Efficient.000000349.081209.CoeffsOnlyFLOAT6040
18 */
19
20 c2 = 0.03462;
21 c3 = 0.09702435063878044594828361677100617517633;
22 c4 = 0.1455365259581706689224254251565092627645;
23 c5 = 0.561;
24 c6 = 0.2290079115904850126662751771814700052182;
25 c7 = 0.5449920884095149873337248228185299947818;
26 c8 = 0.645;
27 c9 = 0.48375;
28 c10 = 0.06757;
29 c11 = 0.2500;
30 c12 = 0.6590650618730998549405331618649220295334;
31 c13 = 0.8206;
32 c14 = 0.9012;
33 c15 = 1.0;
34 c16 = 1.0;
35
36 a21 = c2;
37
38 a31 = -0.0389335438857287327017042687229284478532;
39 a32 = 0.1359578945245091786499878854939346230295;
40
41 a41 = 0.03638413148954266723060635628912731569111;
42 a43 = 0.1091523944686280016918190688673819470733;
43
44 a51 = 2.025763914393969636805657604282571047511;
45 a53 = -7.638023836496292020387602153091964592952;
46 a54 = 6.173259922102322383581944548809393545442;
47
48 a61 = 0.05112275589406060872792270881648288397197;
49 a64 = 0.1770823794555021537929910813839068684087;
50 a65 = 0.00080277624092225014536138698108025283759;
51
52 a71 = 0.1316006357975216279279871693164256985334;
53 a74 = -0.2957276252669636417685183174672273730699;
54 a75 = 0.0878137803564295237421124704053886667082;
55 a76 = 0.62130529752252747743214350056394300261;
56
57 a81 = 0.07166666666666666666666666666666666666667;
58 a86 = 0.3305533578915319409260346730051472207728;
59 a87 = 0.2427799754418013924072986603281861125606;
60
61 a91 = 0.071806640625;
62 a96 = 0.3294380283228177160744825466257672816401;
63 a97 = 0.1165190029271822839255174533742327183599;
64 a98 = -0.034013671875;
65
66 a101 = 0.04836757646340646986611287718844085773549;
67 a106 = 0.03928989925676163974333190042057047002852;
68 a107 = 0.1054740945890344608263649267140088017604;
69 a108 = -0.02143865284648312665982642293830533996214;
70 a109 = -0.1041229174627194437759832813847147895623;
71
72 a111 = -0.02664561487201478635337289243849737340534;
73 a116 = 0.03333333333333333333333333333333333333333;
74 a117 = -0.1631072244872467239162704487554706387141;
75 a118 = 0.03396081684127761199487954930015522928244;
76 a119 = 0.1572319413814626097110769806810024118077;
77 a1110 = 0.2152267478031879552303534778794770376960;
78
79 a121 = 0.03689009248708622334786359863227633989718;
80 a126 = -0.1465181576725542928653609891758501156785;
81 a127 = 0.2242577768172024345345469822625833796001;
82 a128 = 0.02294405717066072637090897902753790803034;
83 a129 = -0.0035850052905728761357394424889330334334;
84 a1210 = 0.08669223316444385506869203619044453906053;
85 a1211 = 0.4383840651968337846196219974168630120572;
86
87 a131 = -0.4866012215113340846662212357570395295088;
88 a136 = -6.304602650282852990657772792012007122988;
89 a137 = -0.281245618289472564778284183790118418111;
90 a138 = -2.679019236219849057687906597489223155566;
91 a139 = 0.518815663924157511565311164615012522024;
92 a1310 = 1.365353187603341710683633635235238678626;
93 a1311 = 5.885091088503946585721274891680604830712;
94 a1312 = 2.802808786272062889819965117517532194812;
95
96 a141 = 0.4185367457753471441471025246471931649633;
97 a146 = 6.724547581906459363100870806514855026676;
98 a147 = -0.425444280164611790606983409697113064616;
99 a148 = 3.343279153001265577811816947557982637749;
100 a149 = 0.617081663117537759528421117507709784737;
101 a1410 = -0.929966123939932833937749523988800852013;
102 a1411 = -6.099948804751010722472962837945508844846;
103 a1412 = -3.002206187889399044804158084895173690015;
104 a1413 = 0.2553202529443445472336424602988558373637;
105
106 a151 = -0.779374086122884664644623040843840506343;
107 a156 = -13.93734253810777678786523664804936051203;
108 a157 = 1.252048853379357320949735183924200895136;
109 a158 = -14.69150040801686878191527989293072091588;
110 a159 = -0.494705058533141685655191992136962873577;
111 a1510 = 2.242974909146236657906984549543692874755;
112 a1511 = 13.36789380382864375813864978592679139881;
113 a1512 = 14.39665048665068644512236935340272139005;
114 a1513 = -0.7975813331776800379127866056663258667437;
115 a1514 = 0.4409353709534277758753793068298041158235;
116
117 a161 = 2.058051337466886442151242368989994043993;
118 a166 = 22.35793772796803295519317565842520212899;
119 a167 = 0.90949810997556332745009198137971890783;
120 a168 = 35.89110098240264104710550686568482456493;
121 a169 = -3.442515027624453437985000403608480262211;
122 a1610 = -4.865481358036368826566013387928704014496;
123 a1611 = -18.90980381354342625688427480879773032857;
124 a1612 = -34.26354448030451782929251177395134170515;
125 a1613 = 1.264756521695642578827783499806516664686;
126
127 b1 = 0.01461197685842315252051541915018784713459;
128 b8 = -0.3915211862331339089410228267288242030810;
129 b9 = 0.2310932500289506415909675644868993669908;
130 b10 = 0.1274766769992852382560589467488989175618;
131 b11 = 0.2246434176204157731566981937082069688984;
132 b12 = 0.5684352689748512932705226972873692126743;
133 b13 = 0.05825871557215827200814768021863420902155;
134 b14 = 0.1364317403482215641609022744494239843327;
135 b15 = 0.03057013983082797397721005067920369646664;
136
137 d1 = 0.01996996514886773085518508418098868756464;
138 d8 = 2.191499304949330054530747099310837524864;
139 d9 = 0.08857071848208438030833722031786358862953;
140 d10 = 0.1140560234865965622484956605091432032674;
141 d11 = 0.2533163805345107065564577734569651977347;
142 d12 = -2.056564386240941011158999594595981300493;
143 d13 = 0.3408096799013119935160094894224543812830;
144 d16 = 0.04834231373823958314376726739772871714902;
145}
146
147void OdeVern98::step_ (double dt) {
148
149 unsigned long i;
150
151 //------------------------------------------------------------------
152 ode_fun_(sol_, k_[0]);
153
154 //------------------------------------------------------------------
155 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*a21*k_[0][i];
156 ode_fun_(soltemp_, k_[1]);
157
158 //------------------------------------------------------------------
159 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a31*k_[0][i]
160 + a32*k_[1][i]);
161 ode_fun_(soltemp_, k_[2]);
162
163 //------------------------------------------------------------------
164 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a41*k_[0][i]
165 + a43*k_[2][i]);
166 ode_fun_(soltemp_, k_[3]);
167
168 //------------------------------------------------------------------
169 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a51*k_[0][i]
170 + a53*k_[2][i]
171 + a54*k_[3][i]);
172 ode_fun_(soltemp_, k_[4]);
173
174 //------------------------------------------------------------------
175 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a61*k_[0][i]
176 + a64*k_[3][i]
177 + a65*k_[4][i]);
178 ode_fun_(soltemp_, k_[5]);
179
180 //------------------------------------------------------------------
181 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a71*k_[0][i]
182 + a74*k_[3][i]
183 + a75*k_[4][i]
184 + a76*k_[5][i]);
185 ode_fun_(soltemp_, k_[6]);
186
187 //------------------------------------------------------------------
188 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a81*k_[0][i]
189 + a86*k_[5][i]
190 + a87*k_[6][i]);
191 ode_fun_(soltemp_, k_[7]);
192
193 //------------------------------------------------------------------
194 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a91*k_[0][i]
195 + a96*k_[5][i]
196 + a97*k_[6][i]
197 + a98*k_[7][i]);
198 ode_fun_(soltemp_, k_[8]);
199
200 //------------------------------------------------------------------
201 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a101*k_[0][i]
202 + a106*k_[5][i]
203 + a107*k_[6][i]
204 + a108*k_[7][i]
205 + a109*k_[8][i]);
206 ode_fun_(soltemp_, k_[9]);
207
208 //------------------------------------------------------------------
209 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a111*k_[0][i]
210 + a116*k_[5][i]
211 + a117*k_[6][i]
212 + a118*k_[7][i]
213 + a119*k_[8][i]
214 + a1110*k_[9][i]);
215 ode_fun_(soltemp_, k_[10]);
216
217 //------------------------------------------------------------------
218 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a121*k_[0][i]
219 + a126*k_[5][i]
220 + a127*k_[6][i]
221 + a128*k_[7][i]
222 + a129*k_[8][i]
223 + a1210*k_[9][i]
224 + a1211*k_[10][i]);
225 ode_fun_(soltemp_, k_[11]);
226
227 //------------------------------------------------------------------
228 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a131*k_[0][i]
229 + a136*k_[5][i]
230 + a137*k_[6][i]
231 + a138*k_[7][i]
232 + a139*k_[8][i]
233 + a1310*k_[9][i]
234 + a1311*k_[10][i]
235 + a1312*k_[11][i]);
236 ode_fun_(soltemp_, k_[12]);
237
238 //------------------------------------------------------------------
239 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a141*k_[0][i]
240 + a146*k_[5][i]
241 + a147*k_[6][i]
242 + a148*k_[7][i]
243 + a149*k_[8][i]
244 + a1410*k_[9][i]
245 + a1411*k_[10][i]
246 + a1412*k_[11][i]
247 + a1413*k_[12][i]);
248 ode_fun_(soltemp_, k_[13]);
249
250 //------------------------------------------------------------------
251 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a151*k_[0][i]
252 + a156*k_[5][i]
253 + a157*k_[6][i]
254 + a158*k_[7][i]
255 + a159*k_[8][i]
256 + a1510*k_[9][i]
257 + a1511*k_[10][i]
258 + a1512*k_[11][i]
259 + a1513*k_[12][i]
260 + a1514*k_[13][i]);
261 ode_fun_(soltemp_, k_[14]);
262
263 //------------------------------------------------------------------
264 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a161*k_[0][i]
265 + a166*k_[5][i]
266 + a167*k_[6][i]
267 + a168*k_[7][i]
268 + a169*k_[8][i]
269 + a1610*k_[9][i]
270 + a1611*k_[10][i]
271 + a1612*k_[11][i]
272 + a1613*k_[12][i]);
273 ode_fun_(soltemp_, k_[15]);
274
275 //------------------------------------------------------------------
276 for (i=0; i<neq_; i++) {
277 solemb_[i] = sol_[i] + dt*(d1*k_[0][i]
278 + d8*k_[7][i]
279 + d9*k_[8][i]
280 + d10*k_[9][i]
281 + d11*k_[10][i]
282 + d12*k_[11][i]
283 + d13*k_[12][i]
284 + d16*k_[15][i]);
285 sol_[i] = sol_[i] + dt*(b1*k_[0][i]
286 + b8*k_[7][i]
287 + b9*k_[8][i]
288 + b10*k_[9][i]
289 + b11*k_[10][i]
290 + b12*k_[11][i]
291 + b13*k_[12][i]
292 + b14*k_[13][i]
293 + b15*k_[14][i]);
294 }
295}
296
297} // 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::OdeVern98::OdeVern98
OdeVern98(unsigned long neq)
constructs
Definition ode_vern_98.cc:7
ode_vern_98.h