src/Propagator/ButcherTableau.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Propagator/ButcherTableau.hpp"
0010
0011 #include <cmath>
0012 #include <numeric>
0013 #include <stdexcept>
0014
0015 namespace Acts {
0016
0017 ButcherTableau::ButcherTableau(std::string name, unsigned order,
0018                                unsigned embeddedOrder, std::vector<double> c,
0019                                std::vector<std::vector<double>> a,
0020                                std::vector<double> b,
0021                                std::vector<double> bEmbedded)
0022     : m_name(std::move(name)),
0023       m_order(order),
0024       m_embeddedOrder(embeddedOrder),
0025       m_c(std::move(c)),
0026       m_b(std::move(b)),
0027       m_bEmbedded(std::move(bEmbedded)) {
0028   const std::size_t s = m_c.size();
0029   if (s == 0 || m_b.size() != s || a.size() != s ||
0030       (!m_bEmbedded.empty() && m_bEmbedded.size() != s)) {
0031     throw std::invalid_argument("ButcherTableau " + m_name +
0032                                 ": the sizes do not match");
0033   }
0034   if (m_c[0] != 0.) {
0035     throw std::invalid_argument("ButcherTableau " + m_name +
0036                                 ": the first node must be zero");
0037   }
0038
0039   m_a.assign(s * s, 0.);
0040   for (std::size_t i = 0; i < s; ++i) {
0041     if (a[i].size() != i) {
0042       throw std::invalid_argument("ButcherTableau " + m_name + ": row " +
0043                                   std::to_string(i) +
0044                                   " must have one coefficient per previous "
0045                                   "stage");
0046     }
0047     // The round-off of the row sum grows with the size of the coefficients,
0048     // so the tolerance is relative to the sum of their magnitudes.
0049     const double rowSum = std::accumulate(a[i].begin(), a[i].end(), 0.);
0050     const double rowNorm = std::accumulate(
0051         a[i].begin(), a[i].end(), 0.,
0052         [](double sum, double aij) { return sum + std::abs(aij); });
0053     if (std::abs(rowSum - m_c[i]) > 1e-13 * (1. + rowNorm)) {
0054       throw std::invalid_argument("ButcherTableau " + m_name + ": row " +
0055                                   std::to_string(i) +
0056                                   " does not sum to its node");
0057     }
0058     for (std::size_t j = 0; j < i; ++j) {
0059       m_a[i * s + j] = a[i][j];
0060     }
0061   }
0062
0063   // Exact equality, so that the last stage is exactly at the end
0064   const std::size_t last = s - 1;
0065   m_firstSameAsLast = s > 1 && m_c[last] == 1. && m_b[last] == 0.;
0066   for (std::size_t j = 0; j < last && m_firstSameAsLast; ++j) {
0067     m_firstSameAsLast = m_a[last * s + j] == m_b[j];
0068   }
0069 }
0070
0071 std::shared_ptr<const ButcherTableau> ButcherTableau::classicalRk4() {
0072   static const auto tableau = std::make_shared<const ButcherTableau>(
0073       "ClassicalRk4", 4, 0, std::vector<double>{0., 1. / 2, 1. / 2, 1.},
0074       std::vector<std::vector<double>>{
0075           {}, {1. / 2}, {0., 1. / 2}, {0., 0., 1.}},
0076       std::vector<double>{1. / 6, 1. / 3, 1. / 3, 1. / 6},
0077       std::vector<double>{});
0078   return tableau;
0079 }
0080
0081 std::shared_ptr<const ButcherTableau> ButcherTableau::dormandPrince54() {
0082   static const auto tableau = std::make_shared<const ButcherTableau>(
0083       "DormandPrince54", 5, 4,
0084       std::vector<double>{0., 1. / 5, 3. / 10, 4. / 5, 8. / 9, 1., 1.},
0085       std::vector<std::vector<double>>{
0086           {},
0087           {1. / 5},
0088           {3. / 40, 9. / 40},
0089           {44. / 45, -56. / 15, 32. / 9},
0090           {19372. / 6561, -25360. / 2187, 64448. / 6561, -212. / 729},
0091           {9017. / 3168, -355. / 33, 46732. / 5247, 49. / 176, -5103. / 18656},
0092           {35. / 384, 0., 500. / 1113, 125. / 192, -2187. / 6784, 11. / 84}},
0093       std::vector<double>{35. / 384, 0., 500. / 1113, 125. / 192, -2187. / 6784,
0094                           11. / 84, 0.},
0095       std::vector<double>{5179. / 57600, 0., 7571. / 16695, 393. / 640,
0096                           -92097. / 339200, 187. / 2100, 1. / 40});
0097   return tableau;
0098 }
0099
0100 std::shared_ptr<const ButcherTableau> ButcherTableau::verner98() {
0101   // Coefficients from J. H. Verner, https://www.sfu.ca/~jverner/, file
0102   // RKV98.IIa.Efficient.000000349.081210.CoeffsOnlyRADandFLOATS.
0103   //
0104   // The exact coefficients are rationals and rationals times sqrt(6). Where
0105   // the numerator or the denominator does not fit into 53 bits, a double
0106   // literal cannot hold the rational exactly, so the code gives the decimal
0107   // value of the rational instead.
0108   const double sqrt6 = std::sqrt(6.);
0109   static const auto tableau = std::make_shared<const ButcherTableau>(
0110       "Verner98", 9, 8,
0111       std::vector<double>{
0112           0., 1731. / 50000, 7630049. / 53810000 - 983539. / 53810000 * sqrt6,
0113           22890147. / 107620000 - 2950617. / 107620000 * sqrt6, 561. / 1000,
0114           387. / 1000 - 129. / 2000 * sqrt6, 387. / 1000 + 129. / 2000 * sqrt6,
0115           129. / 200, 387. / 800, 6757. / 100000, 1. / 4, 0.6590650618730999,
0116           4103. / 5000, 2253. / 2500, 1., 1.},
0117       std::vector<std::vector<double>>{
0118           {},
0119           {1731. / 50000},
0120           {-177968356965557. / 1002427673820000 +
0121                14180534491313. / 250606918455000 * sqrt6,
0122            64021741529527. / 200485534764000 -
0123                7504450763411. / 100242767382000 * sqrt6},
0124           {22890147. / 430480000 - 2950617. / 430480000 * sqrt6, 0.,
0125            68670441. / 430480000 - 8851851. / 430480000 * sqrt6},
0126           {1.1541094857182945 + 0.35585143038209954 * sqrt6, 0.,
0127            -4.244643541842229 - 1.3853417041862812 * sqrt6,
0128            3.651534056123934 + 1.0294902738041818 * sqrt6},
0129           {11380823631. / 157617812000 - 339148869. / 39404453000 * sqrt6, 0.,
0130            0., 0.27509409584108324 - 0.04001311565984247 * sqrt6,
0131            165912282616977. / 4179075230308000 -
0132                33181894472511. / 2089537615154000 * sqrt6},
0133           {26523528363. / 231790900000 + 863255358. / 123138915625 * sqrt6, 0.,
0134            0., -0.04535066278469676 - 0.10221596690492023 * sqrt6,
0135            0.3181366802236169 - 0.09402893012545875 * sqrt6,
0136            -362925891. / 1690350537500 + 857800423623. / 3380701075000 * sqrt6},
0137           {43. / 600, 0., 0., 0., 0., 43. / 150 + 43. / 2400 * sqrt6,
0138            43. / 150 - 43. / 2400 * sqrt6},
0139           {7353. / 102400, 0., 0., 0., 0.,
0140            22833. / 102400 + 8901. / 204800 * sqrt6,
0141            22833. / 102400 - 8901. / 204800 * sqrt6, -3483. / 102400},
0142           {0.04836757646340647, 0., 0., 0., 0.,
0143            0.07238199692289805 - 0.01350979230006338 * sqrt6,
0144            0.07238199692289805 + 0.01350979230006338 * sqrt6,
0145            -0.021438652846483126, -0.10412291746271944},
0146           {-426968570497. / 54394415898750 -
0147                92754382349. / 12087647977500 * sqrt6,
0148            0., 0., 0., 0., 1. / 30,
0149            -0.024935138638394777 - 0.05640851783761993 * sqrt6,
0150            4389715333607. / 309890657317500 +
0151                92754382349. / 11477431752500 * sqrt6,
0152            4990058173976. / 83757096376875 +
0153                371017529396. / 9306344041875 * sqrt6,
0154            0.1757081936006537 + 0.016133382194788114 * sqrt6},
0155           {0.0013154210043370917 + 0.014523298816645871 * sqrt6, 0., 0., 0., 0.,
0156            -0.1470277986219964 + 0.00020806004636011267 * sqrt6,
0157            -0.03674538349491771 + 0.10655409400308863 * sqrt6,
0158            0.06041011360833149 - 0.0152954534911017 * sqrt6,
0159            0.18124121813885546 - 0.0754549897479561 * sqrt6,
0160            0.16148742604165614 - 0.030535009627036813 * sqrt6,
0161            0.43838406519683376},
0162           {-0.9529699587733558 + 0.19039423971301084 * sqrt6, 0., 0., 0., 0.,
0163            -5.729962318480114 - 0.23459593309005522 * sqrt6,
0164            -4.284206994781277 + 1.6342021387456511 * sqrt6,
0165            -2.187855250643843 - 0.200516857448822 * sqrt6,
0166            2.9418084517982352 - 0.9891826634558618 * sqrt6,
0167            2.345886196104345 - 0.40030092446392296 * sqrt6, 5.8850910885039465,
0168            2.8028087862720628},
0169           {1.0566810851455466 - 0.26052133561707513 * sqrt6, 0., 0., 0., 0.,
0170            6.4842136850429215 + 0.09811590253505763 * sqrt6,
0171            4.505949487169823 - 2.013233075118432 * sqrt6,
0172            2.671206837567015 + 0.27437237384412294 * sqrt6,
0173            -2.698361693987128 + 1.3535240826677508 * sqrt6,
0174            -2.271654661242113 + 0.547742051688576 * sqrt6, -6.099948804751011,
0175            -3.002206187889399, 0.2553202529443446},
0176           {-2.215207438880872 + 0.5861765116544452 * sqrt6, 0., 0., 0., 0.,
0177            -13.845681302507453 - 0.03742054273563755 * sqrt6,
0178            -9.394560230624968 + 4.346459957781801 * sqrt6,
0179            -13.17932865245738 - 0.6173415340948993 * sqrt6,
0180            6.965087119532673 - 3.045447403911066 * sqrt6,
0181            5.261792176682924 - 1.2324269886946428 * sqrt6, 13.367893803828643,
0182            14.396650486650687, -0.79758133317768, 0.4409353709534278},
0183           {3.7340160511235347 - 0.6842097292280843 * sqrt6, 0., 0., 0., 0.,
0184            20.39157340157072 + 0.8027648746808281 * sqrt6,
0185            15.196038809510034 - 5.832455817227349 * sqrt6,
0186            34.126030893706144 + 0.7205868462592434 * sqrt6,
0187            -12.149896310859955 + 3.554773523297931 * sqrt6,
0188            -8.389171072898177 + 1.4385403022174306 * sqrt6, -18.909803813543427,
0189            -34.26354448030452, 1.2647565216956427, 0., 0.}},
0190       std::vector<double>{0.014611976858423152, 0., 0., 0., 0., 0., 0.,
0191                           -0.3915211862331339, 0.23109325002895065,
0192                           0.12747667699928525, 0.2246434176204158,
0193                           0.5684352689748513, 0.058258715572158275,
0194                           0.13643174034822156, 0.030570139830827976, 0.},
0195       std::vector<double>{
0196           0.01996996514886773, 0., 0., 0., 0., 0., 0., 2.19149930494933,
0197           0.08857071848208438, 0.11405602348659656, 0.2533163805345107,
0198           -2.056564386240941, 0.340809679901312, 0., 0., 0.048342313738239585});
0199   return tableau;
0200 }
0201
0202 }  // namespace Acts