diff --git a/src/odeint/butcher_tableau.hpp b/src/odeint/butcher_tableau.hpp index 375e460..8d8eebb 100644 --- a/src/odeint/butcher_tableau.hpp +++ b/src/odeint/butcher_tableau.hpp @@ -1,79 +1,84 @@ /*------------------------------------------------------- - Module : odeint - File : butcher_tableau.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_ODEINT_BUTCHERTABLEAU_HPP #define SPECMICP_ODEINT_BUTCHERTABLEAU_HPP //! \file butcher_tableau.hpp Butcher tableau for embedded runge kutta method #include namespace specmicp { namespace odeint { template class ButcherTableau { public: + static const int RK_order = order; + static const int RK_evaluations = s+1; + ButcherTableau(const std::array& aa, const std::array, s>& ba, const std::array& ca, const std::array& csa ): m_a(aa), m_b(ba), m_c(ca), m_cs(csa) {} ButcherTableau(std::array&& aa, std::array, s>&& ba, std::array&& ca, std::array&& csa ): m_a(aa), m_b(ba), m_c(ca), m_cs(csa) {} double a(int i) const {return m_a[i-2];} double b(int i, int j) const {return m_b[i-2][j-1];} double c(int i) const {return m_c[i-1];} double cs(int i) const {return m_cs[i-1];} + private: std::array m_a; std::array, s> m_b; std::array m_c; std::array m_cs; }; const ButcherTableau<5, 5> butcher_cash_karp45({1.0/5.0, 3.0/10.0, 3.0/5.0, 1.0, 7.0/8.0}, { {{1.0/5.0, 0, 0, 0, 0}, {3.0/40.0, 9.0/40.0, 0, 0, 0}, {3.0/10.0, -9.0/10.0, 6.0/5.0, 0, 0}, {-11.0/54.0, 5.0/2.0, -70.0/27.0, 35.0/27.0, 0}, {1631.0/55296.0, 175.0/512.0, 575.0/13824.0, 44275.0/110592.0, 253.0/4096.0} }}, {37.0/378.0, 0.0, 250.0/621.0, 125.0/594.0, 0.0, 512.0/1771.0}, {2825.0/27648.0, 0.0, 18575.0/48384.0, 13525.0/55296.0, 277.0/14336.0, 1.0/4.0} ); - +using ButcherTableauCashKarp_t = ButcherTableau<5, 5>; const ButcherTableau<5, 6> butcher_dormand_prince45( {1.0/5.0, 3.0/10.0, 4.0/5.0, 8.0/9.0, 1.0, 1.0}, { {{1.0/5.0, 0, 0, 0, 0, 0}, {3.0/40.0, 9.0/40.0, 0, 0, 0, 0}, {44.0/45.0, -56.0/15.0, 32.0/9.0, 0, 0, 0}, {19372.0/6561.0, -25360.0/2187.0, 64448.0/6561.0, -212.0/729.0, 0, 0}, {9017.0/3168.0, -355.0/33.0, 46732.0/5247.0, 49.0/176.0, -5103.0/18656.0, 0}, {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0} }}, {5179.0/57600.0, 0.0, 7571.0/16695.0, 393.0/640.0, -92097.0/339200.0, 187.0/2100.0, 1.0/40.0}, {35.0/384.0, 0.0, 500.0/1113.0, 125.0/192.0, -2187.0/6784.0, 11.0/84.0, 0.0} ); +using ButcherTableauDormandPrince_t = ButcherTableau<5, 6>; } // end namespace odeint } // end namespace specmicp #endif // SPECMICP_ODEINT_BUTCHERTABLEAU_HPP diff --git a/src/odeint/runge_kutta_step.hpp b/src/odeint/runge_kutta_step.hpp index e3a41d0..0a48c1c 100644 --- a/src/odeint/runge_kutta_step.hpp +++ b/src/odeint/runge_kutta_step.hpp @@ -1,83 +1,90 @@ /*------------------------------------------------------- - Module : odeint - File : runge_kutta_step - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #ifndef SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP #define SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP //! \file runge_kutta_step algorithm for embedded runge kutta #include "odeint_types.hpp" #include "butcher_tableau.hpp" namespace specmicp { namespace odeint { //! \brief Base class for computing a step of the Range Kutta algorithm -template +template class EmbeddedRungeKuttaStep { public: using vec_t = vector_t; EmbeddedRungeKuttaStep(const rhs_f& rhs, - const ButcherTableau& butcher): + const butcher_tableau_t& butcher): m_rhs(rhs), m_butcher(butcher) {} //! \brief compute a step //! //! \param y initial value //! \param dydx initial value of the right hand site //! \param x value of the integration variable //! \param timestep the integration step //! \param y_out value at the end of the step (y+timestep) //! \param y_err estimation of the truncation error void run_step(const vec_t& y, const vec_t& dydx, scalar_t x, scalar_t timestep, vec_t& y_out, vec_t& y_err); protected: //! Return the value of the right hand side at x for value y void get_rhs(scalar_t x, const vector_t& y, vec_t& dydx) {return m_rhs(x, y, dydx);} //! Return the Butcher tableau - ButcherTableau& butcher() {return m_butcher;} + butcher_tableau_t& butcher() {return m_butcher;} private: rhs_f m_rhs; - ButcherTableau m_butcher; + butcher_tableau_t m_butcher; }; //! Embedded Runge Kutta of fifth order -template -class EmbeddedRungeKuttaStep45: EmbeddedRungeKuttaStep +template +class EmbeddedRungeKuttaStep45: EmbeddedRungeKuttaStep { public: - using base = EmbeddedRungeKuttaStep; + using base = EmbeddedRungeKuttaStep; using vec_t = typename base::vec_t; using base::butcher; using base::get_rhs; - static_assert(s ==5 || s==6, "Template parameter s must be 5 or 6"); + static_assert(butcher_tableau_t::RK_order == 5, "The butcher tableau must correspond to a Runge-Kutta method of order 5"); + static_assert(butcher_tableau_t::RK_evaluations == 6 || butcher_tableau_t::RK_evaluations == 7, "Lol - try again"); EmbeddedRungeKuttaStep45(const rhs_f& rhs, - const ButcherTableau<5, s>& butcher): + const butcher_tableau_t& butcher): base(rhs, butcher) { } void run_step(const vec_t& y, const vec_t& dydx, scalar_t x, scalar_t timestep, vec_t& y_out, vec_t& y_err); }; +template +using DormandPrinceStep = EmbeddedRungeKuttaStep45; + +template +using CashKarpStep = EmbeddedRungeKuttaStep45; + } // end namespace odeint } // end namespace specmicp #include "runge_kutta_step.inl" #endif // SPECMICP_ODEINT_RUNGEKUTTASTEP_HPP diff --git a/src/odeint/runge_kutta_step.inl b/src/odeint/runge_kutta_step.inl index 1d8d335..7df3c24 100644 --- a/src/odeint/runge_kutta_step.inl +++ b/src/odeint/runge_kutta_step.inl @@ -1,69 +1,69 @@ /*------------------------------------------------------- - Module : odeint - File : runge_kutta_step.cpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "runge_kutta_step.hpp" // syntaxic coloration only namespace specmicp { namespace odeint { -template -void EmbeddedRungeKuttaStep45::run_step(const vec_t& y, +template +void EmbeddedRungeKuttaStep45::run_step(const vec_t& y, const vec_t& dydx, scalar_t x, scalar_t timestep, vec_t& y_out, vec_t& y_err) { vec_t ak2, ak3, ak4, ak5, ak6; vec_t ytemp; ytemp = y + butcher().b(2, 1)*timestep*dydx; get_rhs(x + butcher().a(2)*timestep, ytemp, ak2); ytemp = y + timestep*(butcher().b(3,1)*dydx + butcher().b(3, 2)*ak2); get_rhs(x + butcher().a(3)*timestep, ytemp, ak3); ytemp = y + timestep*(butcher().b(4,1)*dydx + butcher().b(4, 2)*ak2 + butcher().b(4, 2)*ak3); get_rhs(x + butcher().a(4)*timestep, ytemp, ak4); ytemp = y + timestep*(butcher().b(5,1)*dydx + butcher().b(5, 2)*ak2 + butcher().b(5, 3)*ak3 + butcher().b(5, 4)*ak4); get_rhs(x + butcher().a(5)*timestep, ytemp, ak5); ytemp = y + timestep*(butcher().b(6,1)*dydx + butcher().b(6, 2)*ak2 + butcher().b(6, 3)*ak3 + butcher().b(6, 4)*ak4 + butcher().b(6, 5)*ak5); get_rhs(x + butcher().a(6)*timestep, ytemp, ak6); - if (s == 5) + if (butcher_tableau_t::RK_evaluations == 6) { y_out = y + timestep*(butcher().c(1)*dydx + butcher().c(2)*ak2 + butcher().c(3)*ak3 + butcher().c(4)*ak4 + butcher().c(5)*ak5 + butcher().c(6)*ak6); y_err = timestep * ( (butcher().c(1) - butcher().cs(1))*dydx + (butcher().c(2) - butcher().cs(2))*ak2 + (butcher().c(3) - butcher().cs(3))*ak3 + (butcher().c(4) - butcher().cs(4))*ak4 + (butcher().c(5) - butcher().cs(5))*ak5 + (butcher().c(6) - butcher().cs(6))*ak6 ); } else { vec_t ak7; ytemp = y + timestep*(butcher().b(7,1)*dydx + butcher().b(7, 2)*ak2 + butcher().b(7, 3)*ak3 + butcher().b(7, 4)*ak4 + butcher().b(7, 5)*ak5 + butcher().b(7, 6)*ak6); get_rhs(x + butcher().a(7)*timestep, ytemp, ak7); y_out = y + timestep*(butcher().c(1)*dydx + butcher().c(2)*ak2 + butcher().c(3)*ak3 + butcher().c(4)*ak4 + butcher().c(5)*ak5 + butcher().c(6)*ak6 + butcher().c(7)*ak7); y_err = timestep * ( (butcher().c(1) - butcher().cs(1))*dydx + (butcher().c(2) - butcher().cs(2))*ak2 + (butcher().c(3) - butcher().cs(3))*ak3 + (butcher().c(4) - butcher().cs(4))*ak4 + (butcher().c(5) - butcher().cs(5))*ak5 + (butcher().c(6) - butcher().cs(6))*ak6 + (butcher().c(7) - butcher().cs(7))*ak7 ); } } } // end namespace odeint } // end namespace specmicp diff --git a/tests/odeint/test_embeddedrungekutta.hpp b/tests/odeint/test_embeddedrungekutta.hpp index 6ceb439..8c7d485 100644 --- a/tests/odeint/test_embeddedrungekutta.hpp +++ b/tests/odeint/test_embeddedrungekutta.hpp @@ -1,72 +1,72 @@ /*------------------------------------------------------- - Module : tests/odeint - File : test_embeddedrungekutta.hpp - Author : Fabien Georget Copyright (c) 2014, Fabien Georget, Princeton University ---------------------------------------------------------*/ #include "cxxtest/TestSuite.h" #include "odeint/runge_kutta_step.hpp" using namespace specmicp::odeint; void rhs1(double x, const double& yin, double& yout) { yout = -2*yin; } void rhs2(double x, const vector_t<2>& y, vector_t<2>& yout) { yout(0) = y(1); yout(1) = -y(0); } class TestEmbeddedRungeKutta : public CxxTest::TestSuite { public: void test_simpletest_cashkarp() { - EmbeddedRungeKuttaStep45<1, 5> algo(rhs1, butcher_cash_karp45); + CashKarpStep<1> algo(rhs1, butcher_cash_karp45); double yout, yerr; double tf = 0.01; algo.run_step(1, -2, 0, tf, yout, yerr); TS_ASSERT_DELTA(yout, std::exp(-2*tf), 1e-3); tf = 0.05; algo.run_step(1, -2, 0, tf, yout, yerr); TS_ASSERT_DELTA(yout, std::exp(-2*tf), 1e-2); } void test_simpletest_dormandprince() { - EmbeddedRungeKuttaStep45<1, 6> algo(rhs1, butcher_dormand_prince45); + DormandPrinceStep<1> algo(rhs1, butcher_dormand_prince45); double yout, yerr; double tf = 0.005; algo.run_step(1, -2, 0, tf, yout, yerr); TS_ASSERT_DELTA(yout, std::exp(-2*tf), 1e-3); tf = 0.01; algo.run_step(1, -2, 0, tf, yout, yerr); TS_ASSERT_DELTA(yout, std::exp(-2*tf), 1e-2); } void test_spring_cashkarp() { - EmbeddedRungeKuttaStep45<2, 5> algo(rhs2, butcher_cash_karp45); + CashKarpStep<2> algo(rhs2, butcher_cash_karp45); vector_t<2> yout, yerr; vector_t<2> y0, dydx; y0 << 1, 0; dydx << 0, -1; double tf=0.01; algo.run_step(y0, dydx, 0, tf, yout, yerr); TS_ASSERT_DELTA(yout(0), std::cos(tf), 1e-3); } };