Skip to content

Commit 14afe78

Browse files
WenyinWeiWenyinWei
committed
move coefficients to a standalone dop853_coefficients.hpp
1 parent 2e3f0f6 commit 14afe78

File tree

2 files changed

+214
-109
lines changed

2 files changed

+214
-109
lines changed

‎include/integrators/ode/dop853.hpp‎

Lines changed: 30 additions & 109 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
1-
21
#pragma once
32
#include <core/adaptive_integrator.hpp>
43
#include <core/state_creator.hpp>
4+
#include <integrators/ode/dop853_coefficients.hpp>
55
#include <cmath>
66
#include <stdexcept>
77

@@ -14,6 +14,7 @@ template<system_state S, can_be_time T>
1414
class DOP853DenseOutputHelper {
1515
public:
1616
using value_type = typename DOP853Integrator<S, T>::value_type;
17+
1718
// Dense output for DOP853: ported from Fortran CONTD8
1819
// CON: continuous output coefficients, size 8*nd
1920
// ICOMP: index mapping, size nd
@@ -234,91 +235,11 @@ class DOP853Integrator : public AdaptiveIntegrator<S, T> {
234235
}
235236

236237
private:
237-
// DOP853 Butcher tableau coefficients (from Fortran code)
238-
static constexpr time_type c2 = 0.0526001519587677318785587544488;
239-
static constexpr time_type c3 = 0.0789002279381515978178381316732;
240-
static constexpr time_type c4 = 0.118350341907227396726757197510;
241-
static constexpr time_type c5 = 0.281649658092772603273242802490;
242-
static constexpr time_type c6 = 0.333333333333333333333333333333;
243-
static constexpr time_type c7 = 0.25;
244-
static constexpr time_type c8 = 0.307692307692307692307692307692;
245-
static constexpr time_type c9 = 0.651282051282051282051282051282;
246-
static constexpr time_type c10 = 0.6;
247-
static constexpr time_type c11 = 0.857142857142857142857142857142;
248-
249-
// a_ij coefficients (only nonzero, as in Fortran)
250-
static constexpr time_type a21 = 0.0526001519587677318785587544488;
251-
static constexpr time_type a31 = 0.0197250569845378994544595329183;
252-
static constexpr time_type a32 = 0.0591751709536136983633785987549;
253-
static constexpr time_type a41 = 0.0295875854768068491816892993775;
254-
static constexpr time_type a43 = 0.0887627564304205475450678981324;
255-
static constexpr time_type a51 = 0.241365134159266685502369798665;
256-
static constexpr time_type a53 = -0.884549479328286085344864962717;
257-
static constexpr time_type a54 = 0.924834003261792003115737966543;
258-
static constexpr time_type a61 = 0.037037037037037037037037037037;
259-
static constexpr time_type a64 = 0.170828608729473871279604482173;
260-
static constexpr time_type a65 = 0.125467687566822425016691814123;
261-
static constexpr time_type a71 = 0.037109375;
262-
static constexpr time_type a74 = 0.170252211019544039314978060272;
263-
static constexpr time_type a75 = 0.0602165389804559606850219397283;
264-
static constexpr time_type a76 = -0.017578125;
265-
static constexpr time_type a81 = 0.0370920001185047927108779319836;
266-
static constexpr time_type a84 = 0.170383925712239993810214054705;
267-
static constexpr time_type a85 = 0.107262030446373284651809199168;
268-
static constexpr time_type a86 = -0.0153194377486244017527936158236;
269-
static constexpr time_type a87 = 0.00827378916381402288758473766002;
270-
static constexpr time_type a91 = 0.624110958716075717114429577812;
271-
static constexpr time_type a94 = -3.36089262944694129406857109825;
272-
static constexpr time_type a95 = -0.868219346841726006818189891453;
273-
static constexpr time_type a96 = 27.5920996994467083049415600797;
274-
static constexpr time_type a97 = 20.1540675504778934086186788979;
275-
static constexpr time_type a98 = -43.4898841810699588477366255144;
276-
static constexpr time_type a101 = 0.477662536438264365890433908527;
277-
static constexpr time_type a104 = -2.48811461997166764192642586468;
278-
static constexpr time_type a105 = -0.590290826836842996371446475743;
279-
static constexpr time_type a106 = 21.2300514481811942347288949897;
280-
static constexpr time_type a107 = 15.2792336328824235832596922938;
281-
static constexpr time_type a108 = -33.2882109689848629194453265587;
282-
static constexpr time_type a109 = -0.0203312017085086261358222928593;
283-
static constexpr time_type a111 = -0.93714243008598732571704021658;
284-
static constexpr time_type a114 = 5.18637242884406370830023853209;
285-
static constexpr time_type a115 = 1.09143734899672957818500254654;
286-
static constexpr time_type a116 = -8.14978701074692612513997267357;
287-
static constexpr time_type a117 = -18.5200656599969598641566180701;
288-
static constexpr time_type a118 = 22.7394870993505042818970056734;
289-
static constexpr time_type a119 = 2.49360555267965238987089396762;
290-
static constexpr time_type a1110 = -3.0467644718982195003823669022;
291-
static constexpr time_type a121 = 2.27331014751653820792359768449;
292-
static constexpr time_type a124 = -10.5344954667372501984066689879;
293-
static constexpr time_type a125 = -2.00087205822486249909675718444;
294-
static constexpr time_type a126 = -17.9589318631187989172765950534;
295-
static constexpr time_type a127 = 27.9488845294199600508499808837;
296-
static constexpr time_type a128 = -2.85899827713502369474065508674;
297-
static constexpr time_type a129 = -8.87285693353062954433549289258;
298-
static constexpr time_type a1210 = 12.3605671757943030647266201528;
299-
static constexpr time_type a1211 = 0.643392746015763530355970484046;
300-
301-
// b_i coefficients (8th order solution)
302-
static constexpr time_type b1 = 0.0542937341165687622380535766363;
303-
static constexpr time_type b6 = 4.45031289275240888144113950566;
304-
static constexpr time_type b7 = 1.89151789931450038304281599044;
305-
static constexpr time_type b8 = -5.8012039600105847814672114227;
306-
static constexpr time_type b9 = 0.31116436695781989440891606237;
307-
static constexpr time_type b10 = -0.152160949662516078556178806805;
308-
static constexpr time_type b11 = 0.201365400804030348374776537501;
309-
static constexpr time_type b12 = 0.0447106157277725905176885569043;
310-
311-
// Error coefficients (5th order)
312-
static constexpr time_type er1 = 0.01312004499419488073250102996;
313-
static constexpr time_type er6 = -12.25156446376204440720569753;
314-
static constexpr time_type er7 = -4.957589496572501915214079952;
315-
static constexpr time_type er8 = 16.64377182454986536961530415;
316-
static constexpr time_type er9 = -0.3503288487499736816886487290;
317-
static constexpr time_type er10 = 0.3341791187130174790297318841;
318-
static constexpr time_type er11 = 0.08192320648511571246570742613;
319-
static constexpr time_type er12 = -0.02235530786388629525884427845;
320-
321-
// For brevity, d coefficients for dense output are omitted here
238+
// Helper functions to access coefficients
239+
static constexpr time_type get_c(int i) { return diffeq::integrators::ode::dop853::C<time_type>[i]; }
240+
static constexpr time_type get_a(int i, int j) { return diffeq::integrators::ode::dop853::A<time_type>[i][j]; }
241+
static constexpr time_type get_b(int i) { return diffeq::integrators::ode::dop853::B<time_type>[i]; }
242+
static constexpr time_type get_e5(int i) { return diffeq::integrators::ode::dop853::E5<time_type>[i]; }
322243

323244
void dop853_step(const state_type& y, state_type& y_new, state_type& error, time_type dt) {
324245
// Allocate all needed k vectors
@@ -331,67 +252,67 @@ class DOP853Integrator : public AdaptiveIntegrator<S, T> {
331252

332253
// k2 = f(t + c2*dt, y + dt*a21*k1)
333254
for (std::size_t i = 0; i < y.size(); ++i)
334-
temp[i] = y[i] + dt * a21 * k[0][i];
335-
this->sys_(t + c2 * dt, temp, k[1]);
255+
temp[i] = y[i] + dt * get_a(1, 0) * k[0][i];
256+
this->sys_(t + get_c(1) * dt, temp, k[1]);
336257

337258
// k3 = f(t + c3*dt, y + dt*(a31*k1 + a32*k2))
338259
for (std::size_t i = 0; i < y.size(); ++i)
339-
temp[i] = y[i] + dt * (a31 * k[0][i] + a32 * k[1][i]);
340-
this->sys_(t + c3 * dt, temp, k[2]);
260+
temp[i] = y[i] + dt * (get_a(2, 0) * k[0][i] + get_a(2, 1) * k[1][i]);
261+
this->sys_(t + get_c(2) * dt, temp, k[2]);
341262

342263
// k4 = f(t + c4*dt, y + dt*(a41*k1 + a43*k3))
343264
for (std::size_t i = 0; i < y.size(); ++i)
344-
temp[i] = y[i] + dt * (a41 * k[0][i] + a43 * k[2][i]);
345-
this->sys_(t + c4 * dt, temp, k[3]);
265+
temp[i] = y[i] + dt * (get_a(3, 0) * k[0][i] + get_a(3, 2) * k[2][i]);
266+
this->sys_(t + get_c(3) * dt, temp, k[3]);
346267

347268
// k5 = f(t + c5*dt, y + dt*(a51*k1 + a53*k3 + a54*k4))
348269
for (std::size_t i = 0; i < y.size(); ++i)
349-
temp[i] = y[i] + dt * (a51 * k[0][i] + a53 * k[2][i] + a54 * k[3][i]);
350-
this->sys_(t + c5 * dt, temp, k[4]);
270+
temp[i] = y[i] + dt * (get_a(4, 0) * k[0][i] + get_a(4, 2) * k[2][i] + get_a(4, 3) * k[3][i]);
271+
this->sys_(t + get_c(4) * dt, temp, k[4]);
351272

352273
// k6 = f(t + c6*dt, y + dt*(a61*k1 + a64*k4 + a65*k5))
353274
for (std::size_t i = 0; i < y.size(); ++i)
354-
temp[i] = y[i] + dt * (a61 * k[0][i] + a64 * k[3][i] + a65 * k[4][i]);
355-
this->sys_(t + c6 * dt, temp, k[5]);
275+
temp[i] = y[i] + dt * (get_a(5, 0) * k[0][i] + get_a(5, 3) * k[3][i] + get_a(5, 4) * k[4][i]);
276+
this->sys_(t + get_c(5) * dt, temp, k[5]);
356277

357278
// k7 = f(t + c7*dt, y + dt*(a71*k1 + a74*k4 + a75*k5 + a76*k6))
358279
for (std::size_t i = 0; i < y.size(); ++i)
359-
temp[i] = y[i] + dt * (a71 * k[0][i] + a74 * k[3][i] + a75 * k[4][i] + a76 * k[5][i]);
360-
this->sys_(t + c7 * dt, temp, k[6]);
280+
temp[i] = y[i] + dt * (get_a(6, 0) * k[0][i] + get_a(6, 3) * k[3][i] + get_a(6, 4) * k[4][i] + get_a(6, 5) * k[5][i]);
281+
this->sys_(t + get_c(6) * dt, temp, k[6]);
361282

362283
// k8 = f(t + c8*dt, y + dt*(a81*k1 + a84*k4 + a85*k5 + a86*k6 + a87*k7))
363284
for (std::size_t i = 0; i < y.size(); ++i)
364-
temp[i] = y[i] + dt * (a81 * k[0][i] + a84 * k[3][i] + a85 * k[4][i] + a86 * k[5][i] + a87 * k[6][i]);
365-
this->sys_(t + c8 * dt, temp, k[7]);
285+
temp[i] = y[i] + dt * (get_a(7, 0) * k[0][i] + get_a(7, 3) * k[3][i] + get_a(7, 4) * k[4][i] + get_a(7, 5) * k[5][i] + get_a(7, 6) * k[6][i]);
286+
this->sys_(t + get_c(7) * dt, temp, k[7]);
366287

367288
// k9 = f(t + c9*dt, y + dt*(a91*k1 + a94*k4 + a95*k5 + a96*k6 + a97*k7 + a98*k8))
368289
for (std::size_t i = 0; i < y.size(); ++i)
369-
temp[i] = y[i] + dt * (a91 * k[0][i] + a94 * k[3][i] + a95 * k[4][i] + a96 * k[5][i] + a97 * k[6][i] + a98 * k[7][i]);
370-
this->sys_(t + c9 * dt, temp, k[8]);
290+
temp[i] = y[i] + dt * (get_a(8, 0) * k[0][i] + get_a(8, 3) * k[3][i] + get_a(8, 4) * k[4][i] + get_a(8, 5) * k[5][i] + get_a(8, 6) * k[6][i] + get_a(8, 7) * k[7][i]);
291+
this->sys_(t + get_c(8) * dt, temp, k[8]);
371292

372293
// k10 = f(t + c10*dt, y + dt*(a101*k1 + a104*k4 + a105*k5 + a106*k6 + a107*k7 + a108*k8 + a109*k9))
373294
for (std::size_t i = 0; i < y.size(); ++i)
374-
temp[i] = y[i] + dt * (a101 * k[0][i] + a104 * k[3][i] + a105 * k[4][i] + a106 * k[5][i] + a107 * k[6][i] + a108 * k[7][i] + a109 * k[8][i]);
375-
this->sys_(t + c10 * dt, temp, k[9]);
295+
temp[i] = y[i] + dt * (get_a(9, 0) * k[0][i] + get_a(9, 3) * k[3][i] + get_a(9, 4) * k[4][i] + get_a(9, 5) * k[5][i] + get_a(9, 6) * k[6][i] + get_a(9, 7) * k[7][i] + get_a(9, 8) * k[8][i]);
296+
this->sys_(t + get_c(9) * dt, temp, k[9]);
376297

377298
// k11 = f(t + c11*dt, y + dt*(a111*k1 + a114*k4 + a115*k5 + a116*k6 + a117*k7 + a118*k8 + a119*k9 + a1110*k10))
378299
for (std::size_t i = 0; i < y.size(); ++i)
379-
temp[i] = y[i] + dt * (a111 * k[0][i] + a114 * k[3][i] + a115 * k[4][i] + a116 * k[5][i] + a117 * k[6][i] + a118 * k[7][i] + a119 * k[8][i] + a1110 * k[9][i]);
380-
this->sys_(t + c11 * dt, temp, k[10]);
300+
temp[i] = y[i] + dt * (get_a(10, 0) * k[0][i] + get_a(10, 3) * k[3][i] + get_a(10, 4) * k[4][i] + get_a(10, 5) * k[5][i] + get_a(10, 6) * k[6][i] + get_a(10, 7) * k[7][i] + get_a(10, 8) * k[8][i] + get_a(10, 9) * k[9][i]);
301+
this->sys_(t + get_c(10) * dt, temp, k[10]);
381302

382303
// k12 = f(t + dt, y + dt*(a121*k1 + a124*k4 + a125*k5 + a126*k6 + a127*k7 + a128*k8 + a129*k9 + a1210*k10 + a1211*k11))
383304
for (std::size_t i = 0; i < y.size(); ++i)
384-
temp[i] = y[i] + dt * (a121 * k[0][i] + a124 * k[3][i] + a125 * k[4][i] + a126 * k[5][i] + a127 * k[6][i] + a128 * k[7][i] + a129 * k[8][i] + a1210 * k[9][i] + a1211 * k[10][i]);
305+
temp[i] = y[i] + dt * (get_a(11, 0) * k[0][i] + get_a(11, 3) * k[3][i] + get_a(11, 4) * k[4][i] + get_a(11, 5) * k[5][i] + get_a(11, 6) * k[6][i] + get_a(11, 7) * k[7][i] + get_a(11, 8) * k[8][i] + get_a(11, 9) * k[9][i] + get_a(11, 10) * k[10][i]);
385306
this->sys_(t + dt, temp, k[11]);
386307

387308
// 8th order solution (y_new)
388309
for (std::size_t i = 0; i < y.size(); ++i) {
389-
y_new[i] = y[i] + dt * (b1 * k[0][i] + b6 * k[5][i] + b7 * k[6][i] + b8 * k[7][i] + b9 * k[8][i] + b10 * k[9][i] + b11 * k[10][i] + b12 * k[11][i]);
310+
y_new[i] = y[i] + dt * (get_b(0) * k[0][i] + get_b(5) * k[5][i] + get_b(6) * k[6][i] + get_b(7) * k[7][i] + get_b(8) * k[8][i] + get_b(9) * k[9][i] + get_b(10) * k[10][i] + get_b(11) * k[11][i]);
390311
}
391312

392313
// 5th order error estimate (embedded)
393314
for (std::size_t i = 0; i < y.size(); ++i) {
394-
error[i] = dt * (er1 * k[0][i] + er6 * k[5][i] + er7 * k[6][i] + er8 * k[7][i] + er9 * k[8][i] + er10 * k[9][i] + er11 * k[10][i] + er12 * k[11][i]);
315+
error[i] = dt * (get_e5(0) * k[0][i] + get_e5(5) * k[5][i] + get_e5(6) * k[6][i] + get_e5(7) * k[7][i] + get_e5(8) * k[8][i] + get_e5(9) * k[9][i] + get_e5(10) * k[10][i] + get_e5(11) * k[11][i]);
395316
}
396317
}
397318
};

0 commit comments

Comments
 (0)