From 5febe408c461e2524df34f2688ab0f0812176e4d Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 30 Aug 2023 09:33:11 +0200 Subject: [PATCH 1/9] Add 5th order Dormand Prince integrator --- .../dormand_prince_5.hpp | 49 ++++++ .../dormand_prince_5_implementation.hpp | 155 ++++++++++++++++++ 2 files changed, 204 insertions(+) create mode 100644 src/library/numerical_integration/dormand_prince_5.hpp create mode 100644 src/library/numerical_integration/dormand_prince_5_implementation.hpp diff --git a/src/library/numerical_integration/dormand_prince_5.hpp b/src/library/numerical_integration/dormand_prince_5.hpp new file mode 100644 index 000000000..c759a0564 --- /dev/null +++ b/src/library/numerical_integration/dormand_prince_5.hpp @@ -0,0 +1,49 @@ +/** + * @file dormand_prince_5.hpp + * @brief Class for 5th order Dormand and Prince method + */ + +#ifndef S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_HPP_ +#define S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_HPP_ + +#include "embedded_runge_kutta.hpp" + +namespace libra::numerical_integration { + +/** + * @class DormandPrince5 + * @brief Class for Classical 5th order Dormand and Prince method + */ +template +class DormandPrince5 : public EmbeddedRungeKutta { + public: + /** + * @fn DormandPrince5 + * @brief Constructor + * @param [in] step_width: Step width + */ + DormandPrince5(const double step_width, const InterfaceOde& ode); + /** + * @fn CalcInterpolationState + * @brief Calculate interpolation state + * @param [in] sigma: Sigma value (0 < sigma < 1) for interpolation + * @return : interpolated state x(t0 + sigma * h) + */ + Vector CalcInterpolationState(const double sigma) const override; + + private: + std::vector> coefficients_; + /** + * @fn CalcInterpolationWeights + * @brief Calculate weights for interpolation + * @param [in] sigma: Sigma value (0 < sigma < 1) for interpolation + * @return : weights for interpolation + */ + std::vector CalcInterpolationWeights(const double sigma) const; +}; + +} // namespace libra::numerical_integration + +#include "dormand_prince_5_implementation.hpp" + +#endif // S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_HPP_ diff --git a/src/library/numerical_integration/dormand_prince_5_implementation.hpp b/src/library/numerical_integration/dormand_prince_5_implementation.hpp new file mode 100644 index 000000000..8ab997ff2 --- /dev/null +++ b/src/library/numerical_integration/dormand_prince_5_implementation.hpp @@ -0,0 +1,155 @@ +/** + * @file dormand_prince_5_implementation.hpp + * @brief Implementation of 5th order Dormand and Prince method + * @note Ref: J. R. Dormand and P. J. Prince, "Runge-Kutta Triples", 1986 + * O. Montenbruck and E. Gill, "State interpolation for on-board navigation systems", 2001 + */ + +#ifndef S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_IMPLEMENTATION_HPP_ +#define S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_IMPLEMENTATION_HPP_ + +#include "dormand_prince_5.hpp" + +namespace libra::numerical_integration { + +template +DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde& ode) : EmbeddedRungeKutta(step_width, ode) { + // p=5th/q=6th order 5th order Dormand and Prince (7-stage) + this->number_of_stages_ = 7; + this->approximation_order_ = 5; + this->nodes_.assign(this->number_of_stages_, 0.0); + this->weights_.assign(this->number_of_stages_, 0.0); + this->higher_order_weights_.assign(this->number_of_stages_, 0.0); + this->rk_matrix_.assign(this->number_of_stages_, std::vector(this->number_of_stages_, 0.0)); + + this->nodes_[1] = 1.0 / 5.0; + this->nodes_[2] = 3.0 / 10.0; + this->nodes_[3] = 4.0 / 5.0; + this->nodes_[4] = 8.0 / 9.0; + this->nodes_[5] = 1.0; + this->nodes_[6] = 1.0; + + this->higher_order_weights_[0] = 5179.0 / 57600.0; + this->higher_order_weights_[1] = 0.0; + this->higher_order_weights_[2] = 7571.0 / 16695.0; + this->higher_order_weights_[3] = 393.0 / 640.0; + this->higher_order_weights_[4] = -92097.0 / 339200.0; + this->higher_order_weights_[5] = 187.0 / 2100.0; + this->higher_order_weights_[6] = 1.0 / 40.0; + + this->weights_[0] = 35.0 / 384.0; + this->weights_[1] = 0.0; + this->weights_[2] = 500.0 / 1113.0; + this->weights_[3] = 125.0 / 192.0; + this->weights_[4] = -2187.0 / 6784.0; + this->weights_[5] = 11.0 / 84.0; + this->weights_[6] = 0.0; + + this->rk_matrix_[1][0] = 1.0 / 5.0; + + this->rk_matrix_[2][0] = 3.0 / 40.0; + this->rk_matrix_[2][1] = 9.0 / 40.0; + + this->rk_matrix_[3][0] = 44.0 / 45.0; + this->rk_matrix_[3][1] = -56.0 / 15.0; + this->rk_matrix_[3][2] = 32.0 / 9.0; + + this->rk_matrix_[4][0] = 19372.0 / 6561.0; + this->rk_matrix_[4][1] = -25630.0 / 2187.0; + this->rk_matrix_[4][2] = 64448.0 / 6561.0; + this->rk_matrix_[4][3] = -212.0 / 729.0; + + this->rk_matrix_[5][0] = 9017.0 / 3168.0; + this->rk_matrix_[5][1] = -355.0 / 33.0; + this->rk_matrix_[5][2] = 46732.0 / 5247.0; + this->rk_matrix_[5][3] = 49.0 / 176.0; + this->rk_matrix_[5][4] = -5103.0 / 18656.0; + + this->rk_matrix_[6][0] = 35.0 / 384.0; + this->rk_matrix_[6][1] = 0.0; + this->rk_matrix_[6][2] = 500.0 / 1113.0; + this->rk_matrix_[6][3] = 125.0 / 192.0; + this->rk_matrix_[6][4] = -2187.0 / 6784.0; + this->rk_matrix_[6][5] = 11.0 / 84.0; + + // Interpolation coefficients + libra::Vector<5> coefficients_temp; + coefficients_temp[0] = 11282082432.0; + coefficients_temp[1] = -32272833064.0; + coefficients_temp[2] = 34969693132.0; + coefficients_temp[3] = -13107642775.0; + coefficients_temp[4] = 157015080.0; + coefficients_temp = 1.0 / 11282082432.0 * coefficients_temp; + coefficients_.push_back(coefficients_temp); + + coefficients_temp *= 0.0; + coefficients_.push_back(coefficients_temp); + + coefficients_temp[0] = 0.0; + coefficients_temp[1] = -1323431896.0; + coefficients_temp[2] = 2074956840.0; + coefficients_temp[3] = -914128567.0; + coefficients_temp[4] = 15701508.0; + coefficients_temp = -100.0 / 32700410799.0 * coefficients_temp; + coefficients_.push_back(coefficients_temp); + + coefficients_temp[0] = 0.0; + coefficients_temp[1] = -889289856.0; + coefficients_temp[2] = 2460397220.0; + coefficients_temp[3] = -1518414297.0; + coefficients_temp[4] = 94209048.0; + coefficients_temp = 25.0 / 5641041216.0 * coefficients_temp; + coefficients_.push_back(coefficients_temp); + + coefficients_temp[0] = 0.0; + coefficients_temp[1] = -259006536.0; + coefficients_temp[2] = 687873124.0; + coefficients_temp[3] = -451824525.0; + coefficients_temp[4] = 52338360.0; + coefficients_temp = -2187.0 / 199316789632.0 * coefficients_temp; + coefficients_.push_back(coefficients_temp); + + coefficients_temp[0] = 0.0; + coefficients_temp[1] = -361440756.0; + coefficients_temp[2] = 946554244.0; + coefficients_temp[3] = -661884105.0; + coefficients_temp[4] = 106151040.0; + coefficients_temp = 11.0 / 2467955532.0 * coefficients_temp; + coefficients_.push_back(coefficients_temp); +} + +template +Vector DormandPrince5::CalcInterpolationState(const double sigma) const { + // Calc k7 (slope after state update) + Vector state_7 = + this->previous_state_ + this->step_width_ * (1.0 / 6.0 * this->slope_[0] + 1.0 / 6.0 * this->slope_[4] + 2.0 / 3.0 * this->slope_[5]); + Vector k7 = this->ode_.DerivativeFunction(this->current_independent_variable_, state_7); + + std::vector interpolation_weights = CalcInterpolationWeights(sigma); + + Vector interpolation_state = this->previous_state_; + for (size_t i = 0; i < this->number_of_stages_ - 1; i++) { + interpolation_state = interpolation_state + (sigma * this->step_width_ * interpolation_weights[i]) * this->slope_[i]; + } + interpolation_state = interpolation_state + sigma * this->step_width_ * (interpolation_weights[6] * k7); + return interpolation_state; +} + +template +std::vector DormandPrince5::CalcInterpolationWeights(const double sigma) const { + std::vector interpolation_weights; + interpolation_weights.assign(this->number_of_stages_, 0.0); + + for (size_t stage = 0; stage < this->number_of_stages_ - 1; stage++) { + for (size_t j = 0; j < 5; j++) { + interpolation_weights[stage] += pow(sigma, j) * coefficients_[stage][j]; + } + } + interpolation_weights[this->number_of_stages_ - 1] = + sigma * (1.0 - sigma) * (8293050.0 * pow(sigma, 2.0) - 82437520.0 * sigma + 44764047.0) / 29380423.0; + return interpolation_weights; +} + +} // namespace libra::numerical_integration + +#endif // S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_IMPLEMENTATION_HPP_ From 59fc983580289803e415596c07fb3b2513b1aaa2 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 30 Aug 2023 09:33:28 +0200 Subject: [PATCH 2/9] Add test for DP5 --- .../test_runge_kutta.cpp | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/src/library/numerical_integration/test_runge_kutta.cpp b/src/library/numerical_integration/test_runge_kutta.cpp index aa6796628..58e7f7d7e 100644 --- a/src/library/numerical_integration/test_runge_kutta.cpp +++ b/src/library/numerical_integration/test_runge_kutta.cpp @@ -5,6 +5,7 @@ #include #include "../orbit/kepler_orbit.hpp" +#include "dormand_prince_5.hpp" #include "numerical_integrator_manager.hpp" #include "ode_examples.hpp" #include "runge_kutta_4.hpp" @@ -181,6 +182,56 @@ TEST(NUMERICAL_INTEGRATION, InterpolationQuadraticRkf) { EXPECT_NEAR(estimated_result, state[0], 1e-6); } +/** + * @brief Test for integration with quadratic function with DP5 + */ +TEST(NUMERICAL_INTEGRATION, IntegrateQuadraticDp5) { + double step_width_s = 0.1; + libra::numerical_integration::ExampleQuadraticOde ode; + libra::numerical_integration::DormandPrince5<1> dp5_ode(step_width_s, ode); + + libra::Vector<1> state = dp5_ode.GetState(); + EXPECT_DOUBLE_EQ(0.0, state[0]); + + size_t step_num = 10000; + for (size_t i = 0; i < step_num; i++) { + dp5_ode.Integrate(); + } + state = dp5_ode.GetState(); + double estimated_result = (step_width_s * step_num) * (step_width_s * step_num); + + EXPECT_NEAR(estimated_result, state[0], 1e-6); +} + +/** + * @brief Test for interpolation with quadratic function with DP5 + */ +TEST(NUMERICAL_INTEGRATION, InterpolationQuadraticDp5) { + double step_width_s = 10.0; + libra::numerical_integration::ExampleQuadraticOde ode; + libra::numerical_integration::DormandPrince5<1> dp5_ode(step_width_s, ode); + + libra::Vector<1> state = dp5_ode.GetState(); + EXPECT_DOUBLE_EQ(0.0, state[0]); + + dp5_ode.Integrate(); + + // Final value + state = dp5_ode.GetState(); + double estimated_result = step_width_s * step_width_s; + EXPECT_NEAR(estimated_result, state[0], 1e-8); + // Interpolation value + double sigma = 0.1; + state = dp5_ode.CalcInterpolationState(sigma); + estimated_result = (step_width_s * sigma) * (step_width_s * sigma); + EXPECT_NEAR(estimated_result, state[0], 1e-8); + + sigma = 0.515; + state = dp5_ode.CalcInterpolationState(sigma); + estimated_result = (step_width_s * sigma) * (step_width_s * sigma); + EXPECT_NEAR(estimated_result, state[0], 1e-8); +} + /** * @brief Test for integration with 1D position and velocity function with RK4 */ From a8799953681c1dbbb7fb7beb3551309fe67581d6 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Wed, 30 Aug 2023 14:28:58 +0200 Subject: [PATCH 3/9] Add 2 body test for DP5 and fix the coefficient value --- .../dormand_prince_5.hpp | 2 +- .../dormand_prince_5_implementation.hpp | 33 +++++----- .../test_runge_kutta.cpp | 63 +++++++++++++++++++ 3 files changed, 81 insertions(+), 17 deletions(-) diff --git a/src/library/numerical_integration/dormand_prince_5.hpp b/src/library/numerical_integration/dormand_prince_5.hpp index c759a0564..9a1c676ab 100644 --- a/src/library/numerical_integration/dormand_prince_5.hpp +++ b/src/library/numerical_integration/dormand_prince_5.hpp @@ -32,7 +32,7 @@ class DormandPrince5 : public EmbeddedRungeKutta { Vector CalcInterpolationState(const double sigma) const override; private: - std::vector> coefficients_; + std::vector> coefficients_; //!< Coefficients to calculate interpolation weights /** * @fn CalcInterpolationWeights * @brief Calculate weights for interpolation diff --git a/src/library/numerical_integration/dormand_prince_5_implementation.hpp b/src/library/numerical_integration/dormand_prince_5_implementation.hpp index 8ab997ff2..3bad490a3 100644 --- a/src/library/numerical_integration/dormand_prince_5_implementation.hpp +++ b/src/library/numerical_integration/dormand_prince_5_implementation.hpp @@ -14,7 +14,7 @@ namespace libra::numerical_integration { template DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde& ode) : EmbeddedRungeKutta(step_width, ode) { - // p=5th/q=6th order 5th order Dormand and Prince (7-stage) + // p=5th/q=4th order 5th order Dormand and Prince (7-stage) this->number_of_stages_ = 7; this->approximation_order_ = 5; this->nodes_.assign(this->number_of_stages_, 0.0); @@ -22,6 +22,7 @@ DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde this->higher_order_weights_.assign(this->number_of_stages_, 0.0); this->rk_matrix_.assign(this->number_of_stages_, std::vector(this->number_of_stages_, 0.0)); + this->nodes_[0] = 0.0; this->nodes_[1] = 1.0 / 5.0; this->nodes_[2] = 3.0 / 10.0; this->nodes_[3] = 4.0 / 5.0; @@ -29,21 +30,21 @@ DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde this->nodes_[5] = 1.0; this->nodes_[6] = 1.0; - this->higher_order_weights_[0] = 5179.0 / 57600.0; - this->higher_order_weights_[1] = 0.0; - this->higher_order_weights_[2] = 7571.0 / 16695.0; - this->higher_order_weights_[3] = 393.0 / 640.0; - this->higher_order_weights_[4] = -92097.0 / 339200.0; - this->higher_order_weights_[5] = 187.0 / 2100.0; - this->higher_order_weights_[6] = 1.0 / 40.0; - - this->weights_[0] = 35.0 / 384.0; + this->weights_[0] = 5179.0 / 57600.0; this->weights_[1] = 0.0; - this->weights_[2] = 500.0 / 1113.0; - this->weights_[3] = 125.0 / 192.0; - this->weights_[4] = -2187.0 / 6784.0; - this->weights_[5] = 11.0 / 84.0; - this->weights_[6] = 0.0; + this->weights_[2] = 7571.0 / 16695.0; + this->weights_[3] = 393.0 / 640.0; + this->weights_[4] = -92097.0 / 339200.0; + this->weights_[5] = 187.0 / 2100.0; + this->weights_[6] = 1.0 / 40.0; + + this->higher_order_weights_[0] = 35.0 / 384.0; + this->higher_order_weights_[1] = 0.0; + this->higher_order_weights_[2] = 500.0 / 1113.0; + this->higher_order_weights_[3] = 125.0 / 192.0; + this->higher_order_weights_[4] = -2187.0 / 6784.0; + this->higher_order_weights_[5] = 11.0 / 84.0; + this->higher_order_weights_[6] = 0.0; this->rk_matrix_[1][0] = 1.0 / 5.0; @@ -55,7 +56,7 @@ DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde this->rk_matrix_[3][2] = 32.0 / 9.0; this->rk_matrix_[4][0] = 19372.0 / 6561.0; - this->rk_matrix_[4][1] = -25630.0 / 2187.0; + this->rk_matrix_[4][1] = -25360.0 / 2187.0; this->rk_matrix_[4][2] = 64448.0 / 6561.0; this->rk_matrix_[4][3] = -212.0 / 729.0; diff --git a/src/library/numerical_integration/test_runge_kutta.cpp b/src/library/numerical_integration/test_runge_kutta.cpp index 58e7f7d7e..28b67fb99 100644 --- a/src/library/numerical_integration/test_runge_kutta.cpp +++ b/src/library/numerical_integration/test_runge_kutta.cpp @@ -64,6 +64,27 @@ TEST(NUMERICAL_INTEGRATION, IntegrateLinearRkf) { EXPECT_NEAR(estimated_result, state[0], 1e-6); } +/** + * @brief Test for integration with nominal linear function with DP5 + */ +TEST(NUMERICAL_INTEGRATION, IntegrateLinearDp5) { + double step_width_s = 0.1; + libra::numerical_integration::ExampleLinearOde ode; + libra::numerical_integration::DormandPrince5<1> dp5_ode(step_width_s, ode); + + libra::Vector<1> state = dp5_ode.GetState(); + EXPECT_DOUBLE_EQ(0.0, state[0]); + + size_t step_num = 10000; + for (size_t i = 0; i < step_num; i++) { + dp5_ode.Integrate(); + } + state = dp5_ode.GetState(); + double estimated_result = step_width_s * step_num; + + EXPECT_NEAR(estimated_result, state[0], 1e-6); +} + /** * @brief Test for integration with nominal linear function with NumericalIntegratorManager/RK4 */ @@ -292,6 +313,36 @@ TEST(NUMERICAL_INTEGRATION, Integrate1dPositionVelocityRkf) { EXPECT_NEAR(estimated_result[1], state[1], 1e-6); } +/** + * @brief Test for integration with 1D position and velocity function with DP5 + */ +TEST(NUMERICAL_INTEGRATION, Integrate1dPositionVelocityDp5) { + double step_width_s = 0.1; + libra::numerical_integration::Example1dPositionVelocityOde ode; + libra::numerical_integration::DormandPrince5<2> dp5_ode(step_width_s, ode); + + libra::Vector<2> initial_state(0.0); + initial_state[0] = 0.0; + initial_state[1] = 0.1; + dp5_ode.SetState(0.0, initial_state); + + libra::Vector<2> state = dp5_ode.GetState(); + EXPECT_DOUBLE_EQ(initial_state[0], state[0]); + EXPECT_DOUBLE_EQ(initial_state[1], state[1]); + + size_t step_num = 10000; + for (size_t i = 0; i < step_num; i++) { + dp5_ode.Integrate(); + } + state = dp5_ode.GetState(); + libra::Vector<2> estimated_result(0.0); + estimated_result[0] = (step_width_s * step_num) * initial_state[1] + initial_state[0]; + estimated_result[1] = initial_state[1]; + + EXPECT_NEAR(estimated_result[0], state[0], 1e-6); + EXPECT_NEAR(estimated_result[1], state[1], 1e-6); +} + /** * @brief Accuracy comparison between RK4 and RKF for integration with 2D two body orbit with small eccentricity */ @@ -300,6 +351,7 @@ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitSmallEccentricity) { libra::numerical_integration::Example2dTwoBodyOrbitOde ode; libra::numerical_integration::RungeKutta4<4> rk4_ode(step_width_s, ode); libra::numerical_integration::RungeKuttaFehlberg<4> rkf_ode(step_width_s, ode); + libra::numerical_integration::DormandPrince5<4> dp5_ode(step_width_s, ode); libra::Vector<4> initial_state(0.0); const double eccentricity = 0.1; @@ -309,21 +361,26 @@ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitSmallEccentricity) { initial_state[3] = sqrt((1.0 + eccentricity) / (1.0 - eccentricity)); rk4_ode.SetState(0.0, initial_state); rkf_ode.SetState(0.0, initial_state); + dp5_ode.SetState(0.0, initial_state); libra::Vector<4> state_rk4 = rk4_ode.GetState(); libra::Vector<4> state_rkf = rkf_ode.GetState(); + libra::Vector<4> state_dp5 = dp5_ode.GetState(); for (size_t i = 0; i < 4; i++) { EXPECT_DOUBLE_EQ(initial_state[i], state_rk4[i]); EXPECT_DOUBLE_EQ(initial_state[i], state_rkf[i]); + EXPECT_DOUBLE_EQ(initial_state[i], state_dp5[i]); } size_t step_num = 200; for (size_t i = 0; i < step_num; i++) { rk4_ode.Integrate(); rkf_ode.Integrate(); + dp5_ode.Integrate(); } state_rk4 = rk4_ode.GetState(); state_rkf = rkf_ode.GetState(); + state_dp5 = dp5_ode.GetState(); // Estimation by Kepler Orbit calculation libra::Vector<3> initial_position(0.0); @@ -348,6 +405,12 @@ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitSmallEccentricity) { EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + + error_tolerance = 2e-6; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); } /** From e2fe76ee0e1ee2d040dae0ca3e57a5e74364c1c6 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Mon, 4 Sep 2023 09:24:46 +0200 Subject: [PATCH 4/9] Fix interpolation --- .../dormand_prince_5_implementation.hpp | 12 ++--- .../test_runge_kutta.cpp | 50 ++++++++++++++++--- 2 files changed, 48 insertions(+), 14 deletions(-) diff --git a/src/library/numerical_integration/dormand_prince_5_implementation.hpp b/src/library/numerical_integration/dormand_prince_5_implementation.hpp index 3bad490a3..9e98ff7fc 100644 --- a/src/library/numerical_integration/dormand_prince_5_implementation.hpp +++ b/src/library/numerical_integration/dormand_prince_5_implementation.hpp @@ -1,7 +1,8 @@ /** * @file dormand_prince_5_implementation.hpp * @brief Implementation of 5th order Dormand and Prince method - * @note Ref: J. R. Dormand and P. J. Prince, "Runge-Kutta Triples", 1986 + * @note Ref: J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae", 1980 + * J. R. Dormand and P. J. Prince, "Runge-Kutta Triples", 1986 * O. Montenbruck and E. Gill, "State interpolation for on-board navigation systems", 2001 */ @@ -121,18 +122,13 @@ DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde template Vector DormandPrince5::CalcInterpolationState(const double sigma) const { - // Calc k7 (slope after state update) - Vector state_7 = - this->previous_state_ + this->step_width_ * (1.0 / 6.0 * this->slope_[0] + 1.0 / 6.0 * this->slope_[4] + 2.0 / 3.0 * this->slope_[5]); - Vector k7 = this->ode_.DerivativeFunction(this->current_independent_variable_, state_7); - std::vector interpolation_weights = CalcInterpolationWeights(sigma); Vector interpolation_state = this->previous_state_; - for (size_t i = 0; i < this->number_of_stages_ - 1; i++) { + for (size_t i = 0; i < this->number_of_stages_; i++) { interpolation_state = interpolation_state + (sigma * this->step_width_ * interpolation_weights[i]) * this->slope_[i]; } - interpolation_state = interpolation_state + sigma * this->step_width_ * (interpolation_weights[6] * k7); + return interpolation_state; } diff --git a/src/library/numerical_integration/test_runge_kutta.cpp b/src/library/numerical_integration/test_runge_kutta.cpp index 28b67fb99..50f6ba8b2 100644 --- a/src/library/numerical_integration/test_runge_kutta.cpp +++ b/src/library/numerical_integration/test_runge_kutta.cpp @@ -344,7 +344,7 @@ TEST(NUMERICAL_INTEGRATION, Integrate1dPositionVelocityDp5) { } /** - * @brief Accuracy comparison between RK4 and RKF for integration with 2D two body orbit with small eccentricity + * @brief Accuracy comparison among integrators for integration with 2D two body orbit with small eccentricity */ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitSmallEccentricity) { double step_width_s = 0.1; @@ -414,13 +414,14 @@ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitSmallEccentricity) { } /** - * @brief Accuracy comparison between RK4 and RKF for integration with 2D two body orbit with high eccentricity + * @brief Accuracy comparison among integrators for integration with 2D two body orbit with high eccentricity */ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitLargeEccentricity) { double step_width_s = 0.01; libra::numerical_integration::Example2dTwoBodyOrbitOde ode; libra::numerical_integration::RungeKutta4<4> rk4_ode(step_width_s, ode); libra::numerical_integration::RungeKuttaFehlberg<4> rkf_ode(step_width_s, ode); + libra::numerical_integration::DormandPrince5<4> dp5_ode(step_width_s, ode); libra::Vector<4> initial_state(0.0); const double eccentricity = 0.9; @@ -430,21 +431,26 @@ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitLargeEccentricity) { initial_state[3] = sqrt((1.0 + eccentricity) / (1.0 - eccentricity)); rk4_ode.SetState(0.0, initial_state); rkf_ode.SetState(0.0, initial_state); + dp5_ode.SetState(0.0, initial_state); libra::Vector<4> state_rk4 = rk4_ode.GetState(); libra::Vector<4> state_rkf = rkf_ode.GetState(); + libra::Vector<4> state_dp5 = dp5_ode.GetState(); for (size_t i = 0; i < 4; i++) { EXPECT_DOUBLE_EQ(initial_state[i], state_rk4[i]); EXPECT_DOUBLE_EQ(initial_state[i], state_rkf[i]); + EXPECT_DOUBLE_EQ(initial_state[i], state_dp5[i]); } size_t step_num = 2000; for (size_t i = 0; i < step_num; i++) { rk4_ode.Integrate(); rkf_ode.Integrate(); + dp5_ode.Integrate(); } state_rk4 = rk4_ode.GetState(); state_rkf = rkf_ode.GetState(); + state_dp5 = dp5_ode.GetState(); // Estimation by Kepler Orbit calculation libra::Vector<3> initial_position(0.0); @@ -469,15 +475,22 @@ TEST(NUMERICAL_INTEGRATION, Integrate2dTwoBodyOrbitLargeEccentricity) { EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + + error_tolerance = 5e-3; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); } /** - * @brief Accuracy comparison between RK4 and RKF for integration with 2D two body orbit with small eccentricity + * @brief Interpolation accuracy comparison among integrators for integration with 2D two body orbit with small eccentricity */ TEST(NUMERICAL_INTEGRATION, Interpolation2dTwoBodyOrbitSmallEccentricity) { double step_width_s = 1.0; libra::numerical_integration::Example2dTwoBodyOrbitOde ode; libra::numerical_integration::RungeKuttaFehlberg<4> rkf_ode(step_width_s, ode); + libra::numerical_integration::DormandPrince5<4> dp5_ode(step_width_s, ode); libra::Vector<4> initial_state(0.0); const double eccentricity = 0.1; @@ -486,17 +499,22 @@ TEST(NUMERICAL_INTEGRATION, Interpolation2dTwoBodyOrbitSmallEccentricity) { initial_state[2] = 0.0; initial_state[3] = sqrt((1.0 + eccentricity) / (1.0 - eccentricity)); rkf_ode.SetState(0.0, initial_state); + dp5_ode.SetState(0.0, initial_state); libra::Vector<4> state_rkf = rkf_ode.GetState(); + libra::Vector<4> state_dp5 = dp5_ode.GetState(); for (size_t i = 0; i < 4; i++) { EXPECT_DOUBLE_EQ(initial_state[i], state_rkf[i]); + EXPECT_DOUBLE_EQ(initial_state[i], state_dp5[i]); } size_t step_num = 1; for (size_t i = 0; i < step_num; i++) { rkf_ode.Integrate(); + dp5_ode.Integrate(); } state_rkf = rkf_ode.GetState(); + state_dp5 = dp5_ode.GetState(); // Estimation by Kepler Orbit calculation libra::Vector<3> initial_position(0.0); @@ -511,31 +529,51 @@ TEST(NUMERICAL_INTEGRATION, Interpolation2dTwoBodyOrbitSmallEccentricity) { KeplerOrbit kepler(1.0, oe); kepler.CalcOrbit((double)(step_num * step_width_s) / (24.0 * 60.0 * 60.0)); - double error_tolerance = 1e-1; + double error_tolerance = 5e-2; EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_rkf[0], error_tolerance); EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + error_tolerance = 5e-2; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); + // Interpolation double sigma = 0.11; state_rkf = rkf_ode.CalcInterpolationState(sigma); + state_dp5 = dp5_ode.CalcInterpolationState(sigma); kepler.CalcOrbit((double)(step_num * step_width_s * sigma) / (24.0 * 60.0 * 60.0)); - error_tolerance = 5e-3; + error_tolerance = 2e-3; EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_rkf[0], error_tolerance); EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + error_tolerance = 2e-3; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); + // Interpolation sigma = 0.79; state_rkf = rkf_ode.CalcInterpolationState(sigma); + state_dp5 = dp5_ode.CalcInterpolationState(sigma); kepler.CalcOrbit((double)(step_num * step_width_s * sigma) / (24.0 * 60.0 * 60.0)); - error_tolerance = 5e-2; + error_tolerance = 3e-2; EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_rkf[0], error_tolerance); EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + + error_tolerance = 3e-2; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); } From 7ff9db8af7605dda5dad4b59483841256f957dca Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Mon, 4 Sep 2023 09:35:42 +0200 Subject: [PATCH 5/9] Add interpolation test for large eccentricity --- .../test_runge_kutta.cpp | 95 +++++++++++++++++++ 1 file changed, 95 insertions(+) diff --git a/src/library/numerical_integration/test_runge_kutta.cpp b/src/library/numerical_integration/test_runge_kutta.cpp index 50f6ba8b2..ed1faa5a1 100644 --- a/src/library/numerical_integration/test_runge_kutta.cpp +++ b/src/library/numerical_integration/test_runge_kutta.cpp @@ -577,3 +577,98 @@ TEST(NUMERICAL_INTEGRATION, Interpolation2dTwoBodyOrbitSmallEccentricity) { EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); } + +/** + * @brief Interpolation accuracy comparison among integrators for integration with 2D two body orbit with large eccentricity + */ +TEST(NUMERICAL_INTEGRATION, Interpolation2dTwoBodyOrbitLargeEccentricity) { + double step_width_s = 0.01; + libra::numerical_integration::Example2dTwoBodyOrbitOde ode; + libra::numerical_integration::RungeKuttaFehlberg<4> rkf_ode(step_width_s, ode); + libra::numerical_integration::DormandPrince5<4> dp5_ode(step_width_s, ode); + + libra::Vector<4> initial_state(0.0); + const double eccentricity = 0.8; + initial_state[0] = 1.0 - eccentricity; + initial_state[1] = 0.0; + initial_state[2] = 0.0; + initial_state[3] = sqrt((1.0 + eccentricity) / (1.0 - eccentricity)); + rkf_ode.SetState(0.0, initial_state); + dp5_ode.SetState(0.0, initial_state); + + libra::Vector<4> state_rkf = rkf_ode.GetState(); + libra::Vector<4> state_dp5 = dp5_ode.GetState(); + for (size_t i = 0; i < 4; i++) { + EXPECT_DOUBLE_EQ(initial_state[i], state_rkf[i]); + EXPECT_DOUBLE_EQ(initial_state[i], state_dp5[i]); + } + + size_t step_num = 1; + for (size_t i = 0; i < step_num; i++) { + rkf_ode.Integrate(); + dp5_ode.Integrate(); + } + state_rkf = rkf_ode.GetState(); + state_dp5 = dp5_ode.GetState(); + + // Estimation by Kepler Orbit calculation + libra::Vector<3> initial_position(0.0); + libra::Vector<3> initial_velocity(0.0); + + // Final value + initial_position[0] = initial_state[0]; + initial_position[1] = initial_state[1]; + initial_velocity[0] = initial_state[2]; + initial_velocity[1] = initial_state[3]; + OrbitalElements oe(1.0, 0.0, initial_position, initial_velocity); + KeplerOrbit kepler(1.0, oe); + kepler.CalcOrbit((double)(step_num * step_width_s) / (24.0 * 60.0 * 60.0)); + + double error_tolerance = 1e-5; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_rkf[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + + error_tolerance = 1e-5; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); + + // Interpolation + double sigma = 0.11; + state_rkf = rkf_ode.CalcInterpolationState(sigma); + state_dp5 = dp5_ode.CalcInterpolationState(sigma); + kepler.CalcOrbit((double)(step_num * step_width_s * sigma) / (24.0 * 60.0 * 60.0)); + + error_tolerance = 1e-5; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_rkf[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + + error_tolerance = 1e-5; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); + + // Interpolation + sigma = 0.79; + state_rkf = rkf_ode.CalcInterpolationState(sigma); + state_dp5 = dp5_ode.CalcInterpolationState(sigma); + kepler.CalcOrbit((double)(step_num * step_width_s * sigma) / (24.0 * 60.0 * 60.0)); + + error_tolerance = 1e-5; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_rkf[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_rkf[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_rkf[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_rkf[3], error_tolerance); + + error_tolerance = 1e-5; + EXPECT_NEAR(kepler.GetPosition_i_m()[0], state_dp5[0], error_tolerance); + EXPECT_NEAR(kepler.GetPosition_i_m()[1], state_dp5[1], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[0], state_dp5[2], error_tolerance); + EXPECT_NEAR(kepler.GetVelocity_i_m_s()[1], state_dp5[3], error_tolerance); +} From b7efa7493d880bd255522295354a60f210abf92c Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Mon, 4 Sep 2023 09:56:11 +0200 Subject: [PATCH 6/9] Fix comment --- src/library/numerical_integration/dormand_prince_5.hpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/library/numerical_integration/dormand_prince_5.hpp b/src/library/numerical_integration/dormand_prince_5.hpp index 9a1c676ab..702a8f9ad 100644 --- a/src/library/numerical_integration/dormand_prince_5.hpp +++ b/src/library/numerical_integration/dormand_prince_5.hpp @@ -1,6 +1,9 @@ /** * @file dormand_prince_5.hpp * @brief Class for 5th order Dormand and Prince method + * @note Ref: J. R. Dormand and P. J. Prince, "A family of embedded Runge-Kutta formulae", 1980 + * J. R. Dormand and P. J. Prince, "Runge-Kutta Triples", 1986 + * O. Montenbruck and E. Gill, "State interpolation for on-board navigation systems", 2001 */ #ifndef S2E_LIBRARY_NUMERICAL_INTEGRATION_DORMAND_PRINCE_5_HPP_ @@ -12,7 +15,7 @@ namespace libra::numerical_integration { /** * @class DormandPrince5 - * @brief Class for Classical 5th order Dormand and Prince method + * @brief Class for 5th order Dormand and Prince method */ template class DormandPrince5 : public EmbeddedRungeKutta { @@ -21,6 +24,7 @@ class DormandPrince5 : public EmbeddedRungeKutta { * @fn DormandPrince5 * @brief Constructor * @param [in] step_width: Step width + * @param [in] ode: Ordinary differential equation */ DormandPrince5(const double step_width, const InterfaceOde& ode); /** From 1855ee6aed5e5b9fe805fb17d4f4412ee43d2adb Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Mon, 4 Sep 2023 11:26:53 +0200 Subject: [PATCH 7/9] Add DP5 to integrator manager --- .../numerical_integrator_manager.hpp | 4 +++ .../test_runge_kutta.cpp | 26 +++++++++++++++++++ 2 files changed, 30 insertions(+) diff --git a/src/library/numerical_integration/numerical_integrator_manager.hpp b/src/library/numerical_integration/numerical_integrator_manager.hpp index b329bfe00..299d8c70d 100644 --- a/src/library/numerical_integration/numerical_integrator_manager.hpp +++ b/src/library/numerical_integration/numerical_integrator_manager.hpp @@ -20,6 +20,7 @@ namespace libra::numerical_integration { enum class NumericalIntegrationMethod { kRk4 = 0, //!< 4th order Runge-Kutta kRkf, //!< Runge-Kutta-Fehlberg + kDp5, //!< 5th order Dormand and Prince }; /** @@ -43,6 +44,9 @@ class NumericalIntegratorManager { case NumericalIntegrationMethod::kRkf: integrator_ = std::make_shared>(step_width, ode); break; + case NumericalIntegrationMethod::kDp5: + integrator_ = std::make_shared>(step_width, ode); + break; default: integrator_ = std::make_shared>(step_width, ode); break; diff --git a/src/library/numerical_integration/test_runge_kutta.cpp b/src/library/numerical_integration/test_runge_kutta.cpp index ed1faa5a1..9724d2b04 100644 --- a/src/library/numerical_integration/test_runge_kutta.cpp +++ b/src/library/numerical_integration/test_runge_kutta.cpp @@ -132,6 +132,32 @@ TEST(NUMERICAL_INTEGRATION, IntegrateLinearNumericalIntegratorManagerRkf) { EXPECT_NEAR(estimated_result, state[0], 1e-6); } +/** + * @brief Test for integration with nominal linear function with NumericalIntegratorManager/DP5 + */ +TEST(NUMERICAL_INTEGRATION, IntegrateLinearNumericalIntegratorManagerDp5) { + double step_width_s = 0.1; + libra::numerical_integration::ExampleLinearOde ode; + libra::numerical_integration::NumericalIntegratorManager<1> numerical_integrator(step_width_s, ode, + libra::numerical_integration::NumericalIntegrationMethod::kDp5); + + libra::Vector<1> state = numerical_integrator.GetIntegrator()->GetState(); + EXPECT_DOUBLE_EQ(0.0, state[0]); + + size_t step_num = 10000; + for (size_t i = 0; i < step_num; i++) { + numerical_integrator.GetIntegrator()->Integrate(); + } + state = numerical_integrator.GetIntegrator()->GetState(); + double estimated_result = step_width_s * step_num; + EXPECT_NEAR(estimated_result, state[0], 1e-6); + + double sigma = 0.1; + state = numerical_integrator.GetIntegrator()->CalcInterpolationState(sigma); + estimated_result = step_width_s * (double(step_num) - 1.0 + sigma); + EXPECT_NEAR(estimated_result, state[0], 1e-6); +} + /** * @brief Test for integration with quadratic function with RK4 */ From 4f4f87d15185806a7843535413e589acbc9c47e4 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Mon, 4 Sep 2023 11:29:50 +0200 Subject: [PATCH 8/9] Fix to add include --- .../numerical_integration/numerical_integrator_manager.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/library/numerical_integration/numerical_integrator_manager.hpp b/src/library/numerical_integration/numerical_integrator_manager.hpp index 299d8c70d..93e2c3d51 100644 --- a/src/library/numerical_integration/numerical_integrator_manager.hpp +++ b/src/library/numerical_integration/numerical_integrator_manager.hpp @@ -8,6 +8,7 @@ #include +#include "dormand_prince_5.hpp" #include "runge_kutta_4.hpp" #include "runge_kutta_fehlberg.hpp" From ac134a7a41e8f82a43c13ae6885ce7109886d2e9 Mon Sep 17 00:00:00 2001 From: Satoshi Ikari Date: Mon, 4 Sep 2023 14:29:14 +0200 Subject: [PATCH 9/9] Initialize slope coefficients in constructor for interpolation --- .../numerical_integration/dormand_prince_5_implementation.hpp | 2 ++ src/library/numerical_integration/runge_kutta_4.hpp | 2 ++ .../runge_kutta_fehlberg_implementation.hpp | 2 ++ 3 files changed, 6 insertions(+) diff --git a/src/library/numerical_integration/dormand_prince_5_implementation.hpp b/src/library/numerical_integration/dormand_prince_5_implementation.hpp index 9e98ff7fc..28348cceb 100644 --- a/src/library/numerical_integration/dormand_prince_5_implementation.hpp +++ b/src/library/numerical_integration/dormand_prince_5_implementation.hpp @@ -118,6 +118,8 @@ DormandPrince5::DormandPrince5(const double step_width, const InterfaceOde coefficients_temp[4] = 106151040.0; coefficients_temp = 11.0 / 2467955532.0 * coefficients_temp; coefficients_.push_back(coefficients_temp); + + this->CalcSlope(); } template diff --git a/src/library/numerical_integration/runge_kutta_4.hpp b/src/library/numerical_integration/runge_kutta_4.hpp index 47abe8f68..0ac8402af 100644 --- a/src/library/numerical_integration/runge_kutta_4.hpp +++ b/src/library/numerical_integration/runge_kutta_4.hpp @@ -40,6 +40,8 @@ class RungeKutta4 : public RungeKutta { this->rk_matrix_[1][0] = this->rk_matrix_[2][1] = 0.5; this->rk_matrix_[3][2] = 1.0; + + this->CalcSlope(); } // We did not implement the interpolation for RK4 diff --git a/src/library/numerical_integration/runge_kutta_fehlberg_implementation.hpp b/src/library/numerical_integration/runge_kutta_fehlberg_implementation.hpp index 2fef2fa31..d8dd17724 100644 --- a/src/library/numerical_integration/runge_kutta_fehlberg_implementation.hpp +++ b/src/library/numerical_integration/runge_kutta_fehlberg_implementation.hpp @@ -59,6 +59,8 @@ RungeKuttaFehlberg::RungeKuttaFehlberg(const double step_width, const Interfa this->rk_matrix_[5][2] = -3544.0 / 2565.0; this->rk_matrix_[5][3] = 1859.0 / 4104.0; this->rk_matrix_[5][4] = -11.0 / 40.0; + + this->CalcSlope(); } template