DiffEq - Modern C++ ODE Integration Library 1.0.0
High-performance C++ library for solving ODEs with async signal processing
Loading...
Searching...
No Matches
dop853_coefficients.hpp
1#pragma once
2#include <array>
3#include <cmath>
4
5namespace diffeq::integrators::ode::dop853 {
6
7// Constants from scipy DOP853 implementation
8constexpr int N_STAGES = 12;
9constexpr int N_STAGES_EXTENDED = 16;
10constexpr int INTERPOLATOR_POWER = 7;
11
12// Time step coefficients (c_i)
13template<typename T>
14constexpr std::array<T, N_STAGES_EXTENDED> C = {
15 T(0.0),
16 T(0.526001519587677318785587544488e-01),
17 T(0.789002279381515978178381316732e-01),
18 T(0.118350341907227396726757197510),
19 T(0.281649658092772603273242802490),
20 T(0.333333333333333333333333333333),
21 T(0.25),
22 T(0.307692307692307692307692307692),
23 T(0.651282051282051282051282051282),
24 T(0.6),
25 T(0.857142857142857142857142857142),
26 T(1.0),
27 T(1.0),
28 T(0.1),
29 T(0.2),
30 T(0.777777777777777777777777777778)
31};
32
33// Butcher tableau coefficients (A matrix)
34template<typename T>
35constexpr std::array<std::array<T, N_STAGES_EXTENDED>, N_STAGES_EXTENDED> A = {{
36 // Row 0: all zeros
37 {{T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
38
39 // Row 1
40 {{T(5.26001519587677318785587544488e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
41
42 // Row 2
43 {{T(1.97250569845378994544595329183e-2), T(5.91751709536136983633785987549e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
44
45 // Row 3
46 {{T(2.95875854768068491816892993775e-2), T(0.0), T(8.87627564304205475450678981324e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
47
48 // Row 4
49 {{T(2.41365134159266685502369798665e-1), T(0.0), T(-8.84549479328286085344864962717e-1), T(9.24834003261792003115737966543e-1), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
50
51 // Row 5
52 {{T(3.7037037037037037037037037037e-2), T(0.0), T(0.0), T(1.70828608729473871279604482173e-1), T(1.25467687566822425016691814123e-1), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
53
54 // Row 6
55 {{T(3.7109375e-2), T(0.0), T(0.0), T(1.70252211019544039314978060272e-1), T(6.02165389804559606850219397283e-2), T(-1.7578125e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
56
57 // Row 7
58 {{T(3.70920001185047927108779319836e-2), T(0.0), T(0.0), T(1.70383925712239993810214054705e-1), T(1.07262030446373284651809199168e-1), T(-1.53194377486244017527936158236e-2), T(8.27378916381402288758473766002e-3), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
59
60 // Row 8
61 {{T(6.24110958716075717114429577812e-1), T(0.0), T(0.0), T(-3.36089262944694129406857109825), T(-8.68219346841726006818189891453e-1), T(2.75920996994467083049415600797e1), T(2.01540675504778934086186788979e1), T(-4.34898841810699588477366255144e1), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
62
63 // Row 9
64 {{T(4.77662536438264365890433908527e-1), T(0.0), T(0.0), T(-2.48811461997166764192642586468), T(-5.90290826836842996371446475743e-1), T(2.12300514481811942347288949897e1), T(1.52792336328824235832596922938e1), T(-3.32882109689848629194453265587e1), T(-2.03312017085086261358222928593e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
65
66 // Row 10
67 {{T(-9.3714243008598732571704021658e-1), T(0.0), T(0.0), T(5.18637242884406370830023853209), T(1.09143734899672957818500254654), T(-8.14978701074692612513997267357), T(-1.85200656599969598641566180701e1), T(2.27394870993505042818970056734e1), T(2.49360555267965238987089396762), T(-3.0467644718982195003823669022), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
68
69 // Row 11
70 {{T(2.27331014751653820792359768449), T(0.0), T(0.0), T(-1.05344954667372501984066689879e1), T(-2.00087205822486249909675718444), T(-1.79589318631187989172765950534e1), T(2.79488845294199600508499808837e1), T(-2.85899827713502369474065508674), T(-8.87285693353062954433549289258), T(1.23605671757943030647266201528e1), T(6.43392746015763530355970484046e-1), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0)}},
71
72 // Row 12 (8th order solution coefficients)
73 {{T(5.42937341165687622380535766363e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(4.45031289275240888144113950566), T(1.89151789931450038304281599044), T(-5.8012039600105847814672114227), T(3.1116436695781989440891606237e-1), T(-1.52160949662516078556178806805e-1), T(2.01365400804030348374776537501e-1), T(4.47106157277725905176885569043e-2), T(0.0), T(0.0), T(0.0), T(0.0)}},
74
75 // Row 13 (5th order solution coefficients)
76 {{T(5.61675022830479523392909219681e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(0.0), T(2.53500210216624811088794765333e-1), T(-2.46239037470802489917441475441e-1), T(-1.24191423263816360469010140626e-1), T(1.5329179827876569731206322685e-1), T(8.20105229563468988491666602057e-3), T(7.56789766054569976138603589584e-3), T(-8.298e-3), T(0.0), T(0.0), T(0.0)}},
77
78 // Row 14
79 {{T(3.18346481635021405060768473261e-2), T(0.0), T(0.0), T(0.0), T(0.0), T(2.83009096723667755288322961402e-2), T(5.35419883074385676223797384372e-2), T(-5.49237485713909884646569340306e-2), T(0.0), T(0.0), T(-1.08347328697249322858509316994e-4), T(3.82571090835658412954920192323e-4), T(-3.40465008687404560802977114492e-4), T(1.41312443674632500278074618366e-1), T(0.0), T(0.0)}},
80
81 // Row 15
82 {{T(-4.28896301583791923408573538692e-1), T(0.0), T(0.0), T(0.0), T(0.0), T(-4.69762141536116384314449447206), T(7.68342119606259904184240953878), T(4.06898981839711007970213554331), T(3.56727187455281109270669543021e-1), T(0.0), T(0.0), T(0.0), T(-1.39902416515901462129418009734e-3), T(2.9475147891527723389556272149), T(-9.15095847217987001081870187138), T(0.0)}}
83}};
84
85// 8th order solution coefficients (B = A[12, :])
86template<typename T>
87constexpr std::array<T, N_STAGES> B = {
88 T(5.42937341165687622380535766363e-2), // b1
89 T(0.0), // b2
90 T(0.0), // b3
91 T(0.0), // b4
92 T(0.0), // b5
93 T(4.45031289275240888144113950566), // b6
94 T(1.89151789931450038304281599044), // b7
95 T(-5.8012039600105847814672114227), // b8
96 T(3.1116436695781989440891606237e-1), // b9
97 T(-1.52160949662516078556178806805e-1), // b10
98 T(2.01365400804030348374776537501e-1), // b11
99 T(4.47106157277725905176885569043e-2) // b12
100};
101
102// 5th order error coefficients (E5)
103template<typename T>
104constexpr std::array<T, N_STAGES> E5 = {
105 T(0.1312004499419488073250102996e-1), // er1
106 T(0.0), // er2
107 T(0.0), // er3
108 T(0.0), // er4
109 T(0.0), // er5
110 T(-0.1225156446376204440720569753e+1), // er6
111 T(-0.4957589496572501915214079952), // er7
112 T(0.1664377182454986536961530415e+1), // er8
113 T(-0.3503288487499736816886487290), // er9
114 T(0.3341791187130174790297318841), // er10
115 T(0.8192320648511571246570742613e-1), // er11
116 T(-0.2235530786388629525884427845e-1) // er12
117};
118
119// 3rd order error coefficients (E3)
120template<typename T>
121constexpr std::array<T, N_STAGES> E3 = {
122 T(5.42937341165687622380535766363e-2 - 0.244094488188976377952755905512), // b1 - 0.244094488188976377952755905512
123 T(0.0), // b2
124 T(0.0), // b3
125 T(0.0), // b4
126 T(0.0), // b5
127 T(4.45031289275240888144113950566), // b6
128 T(1.89151789931450038304281599044), // b7
129 T(-5.8012039600105847814672114227), // b8
130 T(3.1116436695781989440891606237e-1 - 0.733846688281611857341361741547), // b9 - 0.733846688281611857341361741547
131 T(-1.52160949662516078556178806805e-1), // b10
132 T(2.01365400804030348374776537501e-1), // b11
133 T(4.47106157277725905176885569043e-2 - 0.220588235294117647058823529412e-1) // b12 - 0.220588235294117647058823529412e-1
134};
135
136// Dense output coefficients (D matrix)
137template<typename T>
138constexpr std::array<std::array<T, N_STAGES_EXTENDED>, INTERPOLATOR_POWER - 3> D = {{
139 // D[0]
140 {{
141 T(-0.84289382761090128651353491142e+1), T(0.0), T(0.0), T(0.0), T(0.0),
142 T(0.56671495351937776962531783590), T(-0.30689499459498916912797304727e+1),
143 T(0.23846676565120698287728149680e+1), T(0.21170345824450282767155149946e+1),
144 T(-0.87139158377797299206789907490), T(0.22404374302607882758541771650e+1),
145 T(0.63157877876946881815570249290), T(-0.88990336451333310820698117400e-1),
146 T(0.18148505520854727256656404962e+2), T(-0.91946323924783554000451984436e+1),
147 T(-0.44360363875948939664310572000e+1)
148 }},
149
150 // D[1]
151 {{
152 T(0.10427508642579134603413151009e+2), T(0.0), T(0.0), T(0.0), T(0.0),
153 T(0.24228349177525818288430175319e+3), T(0.16520045171727028198505394887e+3),
154 T(-0.37454675472269020279518312152e+3), T(-0.22113666853125306036270938578e+2),
155 T(0.77334326684722638389603898808e+1), T(-0.30674084731089398182061213626e+2),
156 T(-0.93321305264302278729567221706e+1), T(0.15697238121770843886131091075e+2),
157 T(-0.31139403219565177677282850411e+2), T(-0.93529243588444783865713862664e+1),
158 T(0.35816841486394083752465898540e+2)
159 }},
160
161 // D[2]
162 {{
163 T(0.19985053242002433820987653617e+2), T(0.0), T(0.0), T(0.0), T(0.0),
164 T(-0.38703730874935176555105901742e+3), T(-0.18917813819516756882830838328e+3),
165 T(0.52780815920542364900561016686e+3), T(-0.11573902539959630126141871134e+2),
166 T(0.68812326946963000169666922661e+1), T(-0.10006050966910838403183860980e+1),
167 T(0.77771377980534432092869265740), T(-0.27782057523535084065932004339e+1),
168 T(-0.60196695231264120758267380846e+2), T(0.84320405506677161018159903784e+2),
169 T(0.11992291136182789328035130030e+2)
170 }},
171
172 // D[3]
173 {{
174 T(-0.25693933462703749003312586129e+2), T(0.0), T(0.0), T(0.0), T(0.0),
175 T(-0.15418974869023643374053993627e+3), T(-0.23152937917604549567536039109e+3),
176 T(0.35763911791061412378285349910e+3), T(0.93405324183624310003907691704e+2),
177 T(-0.37458323136451633156875139351e+2), T(0.10409964950896230045147246184e+3),
178 T(0.29840293426660503123344363579e+2), T(-0.43533456590011143754432175058e+2),
179 T(0.96324553959188282948394950600e+2), T(-0.39177261675615439165231486172e+2),
180 T(-0.14972683625798562581422125276e+3)
181 }}
182}};
183
184} // namespace diffeq::integrators::ode::dop853