libode
Easy to compile, fast ODE integrators as C++ classes
Loading...
Searching...
No Matches
ode_vern_65.cc
Go to the documentation of this file.
1
2
3#include "ode_vern_65.h"
4
5namespace ode {
6
7OdeVern65::OdeVern65 (unsigned long neq) :
8 OdeEmbedded (neq, false, 5),
9 OdeRK (neq, 9),
10 OdeERK (neq) {
11
12 method_ = "Vern65";
13
14 /*
15 coefficients copied from:
16 http://people.math.sfu.ca/~jverner/
17 http://people.math.sfu.ca/~jverner/RKV65.IIIXb.Efficient.00000144617.081204.CoeffsOnlyFLOAT
18 */
19
20 c2 = 0.06;
21 c3 = 0.09593333333333333333333333333333333333333;
22 c4 = 0.1439;
23 c5 = 0.4973;
24 c6 = 0.9725;
25 c7 = 0.9995;
26 c8 = 1.0;
27 c9 = 1.0;
28
29 a21 = c2;
30
31 a31 = 0.01923996296296296296296296296296296296296;
32 a32 = 0.07669337037037037037037037037037037037037;
33
34 a41 = 0.035975;
35 a43 = 0.107925;
36
37 a51 = 1.318683415233148260919747276431735612861;
38 a53 = -5.042058063628562225427761634715637693344;
39 a54 = 4.220674648395413964508014358283902080483;
40
41 a61 = -41.87259166432751461803757780644346812905;
42 a63 = 159.4325621631374917700365669070346830453;
43 a64 = -122.1192135650100309202516203389242140663;
44 a65 = 5.531743066200053768252631238332999150076;
45
46 a71 = -54.43015693531650433250642051294142461271;
47 a73 = 207.0672513650184644273657173866509835987;
48 a74 = -158.6108137845899991828742424365058599469;
49 a75 = 6.991816585950242321992597280791793907096;
50 a76 = -0.01859723106220323397765171799549294623692;
51
52 a81 = -54.66374178728197680241215648050386959351;
53 a83 = 207.9528062553893734515824816699834244238;
54 a84 = -159.2889574744995071508959805871426654216;
55 a85 = 7.018743740796944434698170760964252490817;
56 a86 = -0.01833878590504572306472782005141738268361;
57 a87 = -0.0005119484997882099077875432497245168395840;
58
59 a91 = 0.03438957868357036009278820124728322386520;
60 a94 = 0.2582624555633503404659558098586120858767;
61 a95 = 0.4209371189673537150642551514069801967032;
62 a96 = 4.405396469669310170148836816197095664891;
63 a97 = -176.4831190242986576151740942499002125029;
64 a98 = 172.3641334014150730294022582711902413315;
65
66 b1 = 0.03438957868357036009278820124728322386520;
67 b4 = 0.2582624555633503404659558098586120858767;
68 b5 = 0.4209371189673537150642551514069801967032;
69 b6 = 4.405396469669310170148836816197095664891;
70 b7 = -176.4831190242986576151740942499002125029;
71 b8 = 172.3641334014150730294022582711902413315;
72
73 d1 = 0.04909967648382489730906854927971225836479;
74 d4 = 0.2251112229516524153401395320539875329485;
75 d5 = 0.4694682253029562039431948525047387412553;
76 d6 = 0.8065792249988867707634161808995217981443;
77 d8 = -0.6071194891777959797672951465256217122488;
78 d9 = 0.05686113944047569241147603178766138153594;
79}
80
81void OdeVern65::step_ (double dt) {
82
83 unsigned long i;
84
85 //------------------------------------------------------------------
86 ode_fun_(sol_, k_[0]);
87
88 //------------------------------------------------------------------
89 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*a21*k_[0][i];
90 ode_fun_(soltemp_, k_[1]);
91
92 //------------------------------------------------------------------
93 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a31*k_[0][i]
94 + a32*k_[1][i]);
95 ode_fun_(soltemp_, k_[2]);
96
97 //------------------------------------------------------------------
98 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a41*k_[0][i]
99 + a43*k_[2][i]);
100 ode_fun_(soltemp_, k_[3]);
101
102 //------------------------------------------------------------------
103 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a51*k_[0][i]
104 + a53*k_[2][i]
105 + a54*k_[3][i]);
106 ode_fun_(soltemp_, k_[4]);
107
108 //------------------------------------------------------------------
109 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a61*k_[0][i]
110 + a63*k_[2][i]
111 + a64*k_[3][i]
112 + a65*k_[4][i]);
113 ode_fun_(soltemp_, k_[5]);
114
115 //------------------------------------------------------------------
116 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a71*k_[0][i]
117 + a73*k_[2][i]
118 + a74*k_[3][i]
119 + a75*k_[4][i]
120 + a76*k_[5][i]);
121 ode_fun_(soltemp_, k_[6]);
122
123 //------------------------------------------------------------------
124 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a81*k_[0][i]
125 + a83*k_[2][i]
126 + a84*k_[3][i]
127 + a85*k_[4][i]
128 + a86*k_[5][i]
129 + a87*k_[6][i]);
130 ode_fun_(soltemp_, k_[7]);
131
132 //------------------------------------------------------------------
133 for (i=0; i<neq_; i++) soltemp_[i] = sol_[i] + dt*(a91*k_[0][i]
134 + a94*k_[3][i]
135 + a95*k_[4][i]
136 + a96*k_[5][i]
137 + a97*k_[6][i]
138 + a98*k_[7][i]);
139 ode_fun_(soltemp_, k_[8]);
140
141 //------------------------------------------------------------------
142 for (i=0; i<neq_; i++) {
143 solemb_[i] = sol_[i] + dt*(d1*k_[0][i]
144 + d4*k_[3][i]
145 + d5*k_[4][i]
146 + d6*k_[5][i]
147 + d8*k_[7][i]
148 + d9*k_[8][i]);
149 sol_[i] = sol_[i] + dt*(b1*k_[0][i]
150 + b4*k_[3][i]
151 + b5*k_[4][i]
152 + b6*k_[5][i]
153 + b7*k_[6][i]
154 + b8*k_[7][i]);
155 }
156}
157
158} // 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::OdeVern65::OdeVern65
OdeVern65(unsigned long neq)
constructs
Definition ode_vern_65.cc:7
ode_vern_65.h