Acts/Propagator/ButcherTableau.hpp

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 #pragma once
0010
0011 #include <cstddef>
0012 #include <memory>
0013 #include <string>
0014 #include <vector>
0015
0016 namespace Acts {
0017
0018 /// @brief Coefficients of an explicit Runge-Kutta method
0019 ///
0020 /// The tableau defines the stages
0021 ///
0022 ///   Y_i = y + h * sum_{j < i} a(i, j) * K_j,  K_i = f(Y_i),
0023 ///
0024 /// the solution y_1 = y + h * sum_i b(i) * K_i and optionally an embedded
0025 /// solution with the weights @ref bEmbedded for an error estimate.
0026 ///
0027 /// The equation of motion is a second-order equation for the position. A
0028 /// tableau with the coefficients A and b, applied to the first-order system of
0029 /// the free parameters, is the same as the general Runge-Kutta-Nyström method
0030 /// with the position coefficients A^2 and b A (Hairer, Nørsett, Wanner,
0031 /// Solving Ordinary Differential Equations I, 2nd ed., Springer 1993, Section
0032 /// II.14).
0033 ///
0034 /// @note The Runge-Kutta-Nyström step of @ref EigenStepper and of the ATLAS
0035 ///       stepper is a different fourth-order method than @ref classicalRk4
0036 ///       on the first-order system, so the two do not give the same result.
0037 class ButcherTableau {
0038  public:
0039   /// Construct and validate a tableau.
0040   ///
0041   /// @param name Name of the method
0042   /// @param order Order of the solution
0043   /// @param embeddedOrder Order of the embedded solution, ignored without
0044   ///        embedded weights
0045   /// @param c Stage nodes
0046   /// @param a Stage coefficients, one row per stage, each row with the
0047   ///        coefficients of the previous stages
0048   /// @param b Weights of the solution
0049   /// @param bEmbedded Weights of the embedded solution, or empty
0050   ///
0051   /// @throws std::invalid_argument if the sizes do not match, if the first
0052   ///         node is not zero, or if a row of @p a does not sum to its node
0053   ButcherTableau(std::string name, unsigned order, unsigned embeddedOrder,
0054                  std::vector<double> c, std::vector<std::vector<double>> a,
0055                  std::vector<double> b, std::vector<double> bEmbedded);
0056
0057   /// @return the name of the method
0058   const std::string& name() const { return m_name; }
0059
0060   /// @return the number of stages
0061   std::size_t stages() const { return m_c.size(); }
0062
0063   /// @return the order of the solution
0064   unsigned order() const { return m_order; }
0065
0066   /// @return the order of the embedded solution
0067   unsigned embeddedOrder() const { return m_embeddedOrder; }
0068
0069   /// @return true if the tableau has embedded weights
0070   bool hasEmbedded() const { return !m_bEmbedded.empty(); }
0071
0072   /// @return true if the last stage is at the end of the step (FSAL), so
0073   ///         its field is the field at the end
0074   bool firstSameAsLast() const { return m_firstSameAsLast; }
0075
0076   /// @param i Stage index
0077   /// @return the node of stage @p i
0078   double c(std::size_t i) const { return m_c[i]; }
0079
0080   /// @param i Stage index
0081   /// @param j Index of a previous stage, j < i
0082   /// @return the coefficient of stage @p j in stage @p i
0083   double a(std::size_t i, std::size_t j) const { return m_a[i * stages() + j]; }
0084
0085   /// @param i Stage index
0086   /// @return the weight of stage @p i in the solution
0087   double b(std::size_t i) const { return m_b[i]; }
0088
0089   /// @param i Stage index
0090   /// @return the weight of stage @p i in the embedded solution
0091   double bEmbedded(std::size_t i) const { return m_bEmbedded[i]; }
0092
0093   /// The classical fourth-order method without an embedded solution
0094   /// (W. Kutta, Z. Math. Phys. 46 (1901) 435-453; Hairer, Nørsett, Wanner,
0095   /// Solving Ordinary Differential Equations I, Section II.1)
0096   /// @return the shared tableau
0097   static std::shared_ptr<const ButcherTableau> classicalRk4();
0098
0099   /// The Dormand-Prince 5(4) pair RK5(4)7M (J. R. Dormand, P. J. Prince,
0100   /// A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 6
0101   /// (1980) 19-26, Table 2; Hairer, Nørsett, Wanner, Solving Ordinary
0102   /// Differential Equations I, Section II.5, DOPRI5). The stepper propagates
0103   /// the fifth-order solution and uses the fourth-order solution for the
0104   /// error estimate.
0105   /// @return the shared tableau
0106   static std::shared_ptr<const ButcherTableau> dormandPrince54();
0107
0108   /// The "most efficient" Verner 9(8) pair. The method is described in
0109   /// J. H. Verner, Numerically optimal Runge-Kutta pairs with interpolants,
0110   /// Numerical Algorithms 53 (2010) 383-396, doi:10.1007/s11075-009-9290-3.
0111   /// The coefficients are from J. H. Verner, https://www.sfu.ca/~jverner/,
0112   /// file RKV98.IIa.Efficient.000000349.081210.CoeffsOnlyRADandFLOATS. The
0113   /// stepper propagates the ninth-order solution and uses the eighth-order
0114   /// solution for the error estimate.
0115   /// @return the shared tableau
0116   static std::shared_ptr<const ButcherTableau> verner98();
0117
0118  private:
0119   std::string m_name;
0120   unsigned m_order;
0121   unsigned m_embeddedOrder;
0122   std::vector<double> m_c;
0123   std::vector<double> m_a;
0124   std::vector<double> m_b;
0125   std::vector<double> m_bEmbedded;
0126   bool m_firstSameAsLast = false;
0127 };
0128
0129 }  // namespace Acts