From 1c86d851b8762bbb8e00baa5f196593addba225c Mon Sep 17 00:00:00 2001 From: Rene107 Date: Mon, 21 Mar 2016 17:38:32 +0100 Subject: [PATCH 1/8] Added Runge Kutta 56 variable step size integrator to coefficient list + Unit test --- .../NumericalIntegrators/CMakeLists.txt | 4 + ...unitTestRungeKuttaFehlberg56Integrator.cpp | 258 ++++++++++++++++++ .../rungeKuttaCoefficients.cpp | 8 + .../rungeKuttaCoefficients.h | 1 + 4 files changed, 271 insertions(+) create mode 100644 Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp diff --git a/Tudat/Mathematics/NumericalIntegrators/CMakeLists.txt b/Tudat/Mathematics/NumericalIntegrators/CMakeLists.txt index d95ea43777..1b4dd4ecb9 100755 --- a/Tudat/Mathematics/NumericalIntegrators/CMakeLists.txt +++ b/Tudat/Mathematics/NumericalIntegrators/CMakeLists.txt @@ -82,6 +82,10 @@ add_executable(test_RungeKuttaFehlberg45Integrator "${SRCROOT}${MATHEMATICSDIR}/ setup_custom_test_program(test_RungeKuttaFehlberg45Integrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators") target_link_libraries(test_RungeKuttaFehlberg45Integrator tudat_numerical_integrators tudat_input_output ${Boost_LIBRARIES}) +add_executable(test_RungeKuttaFehlberg56Integrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp") +setup_custom_test_program(test_RungeKuttaFehlberg56Integrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators") +target_link_libraries(test_RungeKuttaFehlberg56Integrator tudat_numerical_integrators tudat_input_output ${Boost_LIBRARIES}) + add_executable(test_RungeKuttaFehlberg78Integrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp") setup_custom_test_program(test_RungeKuttaFehlberg78Integrator "${SRCROOT}${MATHEMATICSDIR}/NumericalIntegrators") target_link_libraries(test_RungeKuttaFehlberg78Integrator tudat_numerical_integrators tudat_input_output ${Boost_LIBRARIES}) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp new file mode 100644 index 0000000000..57acb94030 --- /dev/null +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp @@ -0,0 +1,258 @@ +/* Copyright (c) 2010-2016, Delft University of Technology + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without modification, are + * permitted provided that the following conditions are met: + * - Redistributions of source code must retain the above copyright notice, this list of + * conditions and the following disclaimer. + * - Redistributions in binary form must reproduce the above copyright notice, this list of + * conditions and the following disclaimer in the documentation and/or other materials + * provided with the distribution. + * - Neither the name of the Delft University of Technology nor the names of its contributors + * may be used to endorse or promote products derived from this software without specific + * prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS + * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF + * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE + * COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE + * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED + * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED + * OF THE POSSIBILITY OF SUCH DAMAGE. + * + * Changelog + * YYMMDD Author Comment + * 120203 B. Tong Minh Copied RungeKutta4Stepsize unit test. + * 120207 K. Kumar Adapted to use modified benchmark functions in Tudat Core. + * 120213 K. Kumar Modified getCurrentInterval() to getIndependentVariable(); + * transferred to Boost unit test framework. + * 120321 K. Kumar Updated (Burden and Faires, 2011) benchmark function call. + * 120323 K. Kumar Rewrote unit tests to use benchmark data from + * (Burden and Faires, 2011); removed test against benchmark + * functions; renamed file to RFK45 unit test. Other unit tests + * for other Runge-Kutta methods will appear in other dedicated + * unit test files. + * 120327 K. Kumar Added missing comments; added unit test based on output data + * generated by (The Mathworks, 2012). + * 120328 K. Kumar Moved (Burden and Faires, 2011) test class to its own file; + * modified MATLAB unit tests to test forward and backwards in + * time and "forced" and "free" adaptive step size adjustment and + * moved to separate file (added call to function to run Matlab + * tests); added rollback tests for all cases. + * 120404 K. Kumar Updated Matlab unit test by adding discrete-event data file. + * 130118 K. Kumar Rewrote unit test to make use of testing code for numerical + * integrators migrated to Tudat Core. + * 130906 K. Kumar Updated error tolerances for MuPAD-based tests. + * + * 160321 R. Hoogendoorn Created Runge Kutta 56 test + * + * References + * Burden, R.L., Faires, J.D. Numerical Analysis, 7th Edition, Books/Cole, 2001. + * Montenbruck, O., Gill, E. Satellite Orbits: Models, Methods, Applications, Springer, 2005. + * The MathWorks, Inc. RKF54b, Symbolic Math Toolbox, 2012. + * + * Notes + * For the tests using data from the Symbolic Math Toolbox (MathWorks, 2012), the single step + * and full integration error tolerances were picked to be as small as possible, without + * causing the tests to fail. These values are not deemed to indicate any bugs in the code; + * however, it is important to take these discrepancies into account when using this numerical + * integrator. + * + */ + +#define BOOST_TEST_MAIN + +#include +#include +#include + +#include "Tudat/Mathematics/NumericalIntegrators/rungeKuttaVariableStepSizeIntegrator.h" +#include "Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h" +#include "Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h" + +#include "Tudat/InputOutput/basicInputOutput.h" + +#include +#include + +#include + +namespace tudat +{ +namespace unit_tests +{ + +BOOST_AUTO_TEST_SUITE( test_runge_kutta_fehlberg_56_integrator ) + +using numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative; + +//! Test Compare with Runge Kutta 78 +BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) +{ + std::cout << "Test 1 : Comparison RK56 and RK78" << std::endl; + + // Statederivative function. + using tudat::unit_tests::numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative ; + + // Setup integrator + tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + + tudat::numerical_integrators::RungeKuttaCoefficients Coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + + // Integrator settings + double MinimumStepSize = std::numeric_limits::epsilon() ; + double MaximumStepSize = std::numeric_limits::infinity() ; + double InitialStepSize = 1E-8 ; + double RelativeTolerance = 1E-15 ; + double AbsoluteTolerance = 1E-15 ; + + // Initial conditions + double InitialTime = 0.5 ; + Eigen::VectorXd InitialState(1) ; + InitialState << 1.0 ; + + // Setup integrator + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( + Coeff56, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, + MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( + Coeff78, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, + MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + + double EndTime = 1.5 ; + Eigen::VectorXd Solution56 = integrator56.integrateTo(EndTime,InitialStepSize) ; + Eigen::VectorXd Solution78 = integrator78.integrateTo(EndTime,InitialStepSize) ; + + Eigen::VectorXd Difference = Solution78 - Solution56 ; + + std::cout << "Solution = " << Solution78 << std::endl; + std::cout << "Difference between RK56 and RK78: " << std::endl; + std::cout << Difference << std::endl; + + BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-14 ) ; + + // Error order check + // 78 -> h^7 , 56 -> h^5 +// std::cout << "Next stepsize 56: " << integrator56.getNextStepSize() << std::endl; +// std::cout << "Next stepsize 78: " << integrator78.getNextStepSize() << std::endl; +// std::cout << "h^5 : " << std::pow( integrator56.getNextStepSize() ,5.0) << std::endl; +// std::cout << "h^7 : " << std::pow( integrator78.getNextStepSize() ,7.0) << std::endl; +// std::cout << "h^7 - h^5 : " << std::pow( integrator56.getNextStepSize() ,5.0) - std::pow( integrator78.getNextStepSize() ,7.0) << std::endl; +} + +//! Test Compare with Runge Kutta 78 +BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_v2 ) +{ + std::cout << "" << std::endl; + std::cout << "Test 2 : Comparison RK56 and RK78" << std::endl; + + // Statederivative function. + using tudat::unit_tests::numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative ; + + // Setup integrator + tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + + tudat::numerical_integrators::RungeKuttaCoefficients Coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + + // Integrator settings + double MinimumStepSize = std::numeric_limits::epsilon() ; + double MaximumStepSize = std::numeric_limits::infinity() ; + double InitialStepSize = 1E-8 ; + double RelativeTolerance = 1E-10 ; + double AbsoluteTolerance = 1E-10 ; + + // Initial conditions + double InitialTime = 0.2 ; + Eigen::VectorXd InitialState(1) ; + InitialState << -1.0 ; + + // Setup integrator + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( + Coeff56, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, + MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( + Coeff78, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, + MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + + double EndTime = 2.0 ; + Eigen::VectorXd Solution56 = integrator56.integrateTo(EndTime,InitialStepSize) ; + Eigen::VectorXd Solution78 = integrator78.integrateTo(EndTime,InitialStepSize) ; + + Eigen::VectorXd Difference = Solution78 - Solution56 ; + + std::cout << "Solution = " << Solution78 << std::endl; + std::cout << "Difference between RK56 and RK78: " << std::endl; + std::cout << Difference << std::endl; + + BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-8 ) ; +} + +//! Test Compare with Runge Kutta 78 +BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_VanDerPol ) +{ + std::cout << "" << std::endl; + std::cout << "Test 3 : Comparison RK56 and RK78 (Model: Van Der Pol)" << std::endl; + + // Statederivative function. + using tudat::unit_tests::numerical_integrator_test_functions::computeVanDerPolStateDerivative ; + + // Setup integrator + tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + + tudat::numerical_integrators::RungeKuttaCoefficients Coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + + // Integrator settings + double MinimumStepSize = std::numeric_limits::epsilon() ; + double MaximumStepSize = std::numeric_limits::infinity() ; + double InitialStepSize = 1E-8 ; + double RelativeTolerance = 1E-15 ; + double AbsoluteTolerance = 1E-15 ; + + // Initial conditions + double InitialTime = 0.2 ; + Eigen::VectorXd InitialState(2) ; + InitialState << -1.0 , 1.0 ; + + // Setup integrator + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( + Coeff56, computeVanDerPolStateDerivative, InitialTime , InitialState, MinimumStepSize, + MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( + Coeff78, computeVanDerPolStateDerivative, InitialTime , InitialState, MinimumStepSize, + MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + + double EndTime = 1.4 ; + Eigen::VectorXd Solution56 = integrator56.integrateTo(EndTime,InitialStepSize) ; + Eigen::VectorXd Solution78 = integrator78.integrateTo(EndTime,InitialStepSize) ; + + Eigen::VectorXd Difference = Solution78 - Solution56 ; + + std::cout << "Solution = " << Solution78 << std::endl; + std::cout << "Difference between RK56 and RK78: " << std::endl; + std::cout << Difference << std::endl; + + BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs(Difference(1)) , 1E-13 ) ; +} + +BOOST_AUTO_TEST_SUITE_END( ) + +} // namespace unit_tests +} // namespace tudat diff --git a/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.cpp b/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.cpp index 26a5479dd4..e1090a14da 100644 --- a/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.cpp @@ -441,6 +441,7 @@ const RungeKuttaCoefficients& RungeKuttaCoefficients::get( RungeKuttaCoefficients::CoefficientSets coefficientSet ) { static RungeKuttaCoefficients rungeKuttaFehlberg45Coefficients, + rungeKuttaFehlberg56Coefficients, rungeKuttaFehlberg78Coefficients, rungeKutta87DormandPrinceCoefficients; @@ -453,6 +454,13 @@ const RungeKuttaCoefficients& RungeKuttaCoefficients::get( } return rungeKuttaFehlberg45Coefficients; + case rungeKuttaFehlberg56: + if ( rungeKuttaFehlberg56Coefficients.higherOrder != 6 ) + { + initializeRungeKuttaFehlberg56Coefficients( rungeKuttaFehlberg56Coefficients ); + } + return rungeKuttaFehlberg56Coefficients; + case rungeKuttaFehlberg78: if ( rungeKuttaFehlberg78Coefficients.higherOrder != 8 ) { diff --git a/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h b/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h index 81138605d9..86d1e63f27 100644 --- a/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h +++ b/Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h @@ -119,6 +119,7 @@ struct RungeKuttaCoefficients enum CoefficientSets { rungeKuttaFehlberg45, + rungeKuttaFehlberg56, rungeKuttaFehlberg78, rungeKutta87DormandPrince }; From bbca5e6713af428caf45a2d5b28fb0f130931a4e Mon Sep 17 00:00:00 2001 From: Jacco Geul Date: Tue, 22 Mar 2016 13:06:27 +0100 Subject: [PATCH 2/8] UNFINISHED Add test case from Fehlberg1968. --- ...unitTestRungeKuttaFehlberg56Integrator.cpp | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp index 57acb94030..af71c063ec 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp @@ -68,6 +68,8 @@ #include #include +#include + #include "Tudat/Mathematics/NumericalIntegrators/rungeKuttaVariableStepSizeIntegrator.h" #include "Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h" #include "Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h" @@ -86,8 +88,53 @@ namespace unit_tests BOOST_AUTO_TEST_SUITE( test_runge_kutta_fehlberg_56_integrator ) +Eigen::Vector2d FehlbergLogirithmicTestDifferentialEquation( double x, Eigen::Vector2d state){ + Eigen::Vector2d stateDerivative; + stateDerivative( 0 ) = -2*x*state(0)*log(state(1)); + stateDerivative( 1 ) = 2*x*state(1)*log(state(0)); + return stateDerivative; +} + using numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative; +//! Test Compare with Runge Kutta 78 +BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) +{ + tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + + // Integrator settings + double minimumStepSize = std::numeric_limits::epsilon() ; + double maximumStepSize = std::numeric_limits::infinity() ; + double initialStepSize = 1E-16; + double relativeTolerance = 1E-16; + double absoluteTolerance = 1E-16; + + // Initial conditions + double initialTime = 0.0; + double finalTime = 5.0; + Eigen::Vector2d initialState(exp(1.0), 1.0); + + // Setup integrator + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( + coeff56, boost::bind(FehlbergLogirithmicTestDifferentialEquation,_1,_2), + initialTime , initialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); + + double finalTimeSquared = finalTime * finalTime; + + Eigen::Vector2d numSol = integrator56.integrateTo(finalTime, initialStepSize); + Eigen::Vector2d anaSol(exp(cos(finalTimeSquared)), exp(sin(finalTimeSquared))); + + Eigen::Vector2d computedDifference = numSol - anaSol; + std::cout << computedDifference << std::endl; + // std::cout << solution56 << std::endl; +// BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-14 ) ; + +} + + //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) { From bfb125515ff40f0759d3d3b6181c53f4c9ea2843 Mon Sep 17 00:00:00 2001 From: Rene107 Date: Wed, 23 Mar 2016 08:30:07 +0100 Subject: [PATCH 3/8] Update of unit test for RK56 integrator. Added extra test that compares with results of Fehlberg and with an analytical solution --- ...unitTestRungeKuttaFehlberg56Integrator.cpp | 105 +++++++----------- 1 file changed, 38 insertions(+), 67 deletions(-) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp index af71c063ec..134b0c00eb 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp @@ -64,23 +64,20 @@ #define BOOST_TEST_MAIN -#include -#include #include +#include +#include #include +#include + +#include + #include "Tudat/Mathematics/NumericalIntegrators/rungeKuttaVariableStepSizeIntegrator.h" #include "Tudat/Mathematics/NumericalIntegrators/rungeKuttaCoefficients.h" #include "Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h" -#include "Tudat/InputOutput/basicInputOutput.h" - -#include -#include - -#include - namespace tudat { namespace unit_tests @@ -88,16 +85,12 @@ namespace unit_tests BOOST_AUTO_TEST_SUITE( test_runge_kutta_fehlberg_56_integrator ) -Eigen::Vector2d FehlbergLogirithmicTestDifferentialEquation( double x, Eigen::Vector2d state){ - Eigen::Vector2d stateDerivative; - stateDerivative( 0 ) = -2*x*state(0)*log(state(1)); - stateDerivative( 1 ) = 2*x*state(1)*log(state(0)); - return stateDerivative; -} - using numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative; +using numerical_integrator_test_functions::computeVanDerPolStateDerivative ; +using numerical_integrator_test_functions::computeFehlbergLogirithmicTestODEStateDerivative ; +using numerical_integrator_test_functions::computeAnalyticalStateFehlbergODE ; -//! Test Compare with Runge Kutta 78 +//! Compare with analytical solution BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) { tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = @@ -107,7 +100,7 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) // Integrator settings double minimumStepSize = std::numeric_limits::epsilon() ; double maximumStepSize = std::numeric_limits::infinity() ; - double initialStepSize = 1E-16; + double initialStepSize = 1E-6; // Don't make this too small double relativeTolerance = 1E-16; double absoluteTolerance = 1E-16; @@ -118,31 +111,41 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) // Setup integrator tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - coeff56, boost::bind(FehlbergLogirithmicTestDifferentialEquation,_1,_2), + coeff56, computeFehlbergLogirithmicTestODEStateDerivative , initialTime , initialState, minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - double finalTimeSquared = finalTime * finalTime; + // Obtain numerical solution Eigen::Vector2d numSol = integrator56.integrateTo(finalTime, initialStepSize); - Eigen::Vector2d anaSol(exp(cos(finalTimeSquared)), exp(sin(finalTimeSquared))); - Eigen::Vector2d computedDifference = numSol - anaSol; - std::cout << computedDifference << std::endl; - // std::cout << solution56 << std::endl; -// BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-14 ) ; + // Analytical solution + // (page 30, Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control) + Eigen::Vector2d anaSol = computeAnalyticalStateFehlbergODE(finalTime, initialState) ; + + Eigen::Vector2d computedError = numSol - anaSol; + BOOST_CHECK_SMALL( std::fabs(computedError(0)) , 1E-12 ) ; + BOOST_CHECK_SMALL( std::fabs(computedError(1)) , 1E-12 ) ; + + // Error calculated by -> Fehlberg, E. (1968) page 30 + // Initial stepsize unknown.. + Eigen::VectorXd FehlbergError(2) ; + FehlbergError << 0.1072E-12 , -0.2190E-12 ; + + // sign check -> same sign -> positive + // Not always same sign -> initial step size = 1 or 1E-2, failure computedError(1)/FehlbergError(1) + BOOST_CHECK_GE(computedError(0)/FehlbergError(0) ,0.0); + BOOST_CHECK_GE(computedError(1)/FehlbergError(1) ,0.0); + // Check error is similar + BOOST_CHECK_SMALL(std::fabs(computedError(0)/FehlbergError(0)), 3.0); + BOOST_CHECK_SMALL(std::fabs(computedError(1)/FehlbergError(1)), 1.0); } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) { - std::cout << "Test 1 : Comparison RK56 and RK78" << std::endl; - - // Statederivative function. - using tudat::unit_tests::numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative ; - // Setup integrator tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = tudat::numerical_integrators::RungeKuttaCoefficients::get( @@ -155,14 +158,14 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) // Integrator settings double MinimumStepSize = std::numeric_limits::epsilon() ; double MaximumStepSize = std::numeric_limits::infinity() ; - double InitialStepSize = 1E-8 ; + double InitialStepSize = 1E-4; // Don't make this too small double RelativeTolerance = 1E-15 ; double AbsoluteTolerance = 1E-15 ; // Initial conditions double InitialTime = 0.5 ; Eigen::VectorXd InitialState(1) ; - InitialState << 1.0 ; + InitialState << 0.5 ; // 1 large error // Setup integrator tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( @@ -179,30 +182,12 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) Eigen::VectorXd Difference = Solution78 - Solution56 ; - std::cout << "Solution = " << Solution78 << std::endl; - std::cout << "Difference between RK56 and RK78: " << std::endl; - std::cout << Difference << std::endl; - - BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-14 ) ; - - // Error order check - // 78 -> h^7 , 56 -> h^5 -// std::cout << "Next stepsize 56: " << integrator56.getNextStepSize() << std::endl; -// std::cout << "Next stepsize 78: " << integrator78.getNextStepSize() << std::endl; -// std::cout << "h^5 : " << std::pow( integrator56.getNextStepSize() ,5.0) << std::endl; -// std::cout << "h^7 : " << std::pow( integrator78.getNextStepSize() ,7.0) << std::endl; -// std::cout << "h^7 - h^5 : " << std::pow( integrator56.getNextStepSize() ,5.0) - std::pow( integrator78.getNextStepSize() ,7.0) << std::endl; + BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-13 ) ; } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_v2 ) { - std::cout << "" << std::endl; - std::cout << "Test 2 : Comparison RK56 and RK78" << std::endl; - - // Statederivative function. - using tudat::unit_tests::numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative ; - // Setup integrator tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = tudat::numerical_integrators::RungeKuttaCoefficients::get( @@ -215,7 +200,7 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_v2 ) // Integrator settings double MinimumStepSize = std::numeric_limits::epsilon() ; double MaximumStepSize = std::numeric_limits::infinity() ; - double InitialStepSize = 1E-8 ; + double InitialStepSize = 1; // Don't make this too small double RelativeTolerance = 1E-10 ; double AbsoluteTolerance = 1E-10 ; @@ -239,22 +224,12 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_v2 ) Eigen::VectorXd Difference = Solution78 - Solution56 ; - std::cout << "Solution = " << Solution78 << std::endl; - std::cout << "Difference between RK56 and RK78: " << std::endl; - std::cout << Difference << std::endl; - BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-8 ) ; } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_VanDerPol ) { - std::cout << "" << std::endl; - std::cout << "Test 3 : Comparison RK56 and RK78 (Model: Van Der Pol)" << std::endl; - - // Statederivative function. - using tudat::unit_tests::numerical_integrator_test_functions::computeVanDerPolStateDerivative ; - // Setup integrator tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = tudat::numerical_integrators::RungeKuttaCoefficients::get( @@ -267,7 +242,7 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_VanDerPol ) // Integrator settings double MinimumStepSize = std::numeric_limits::epsilon() ; double MaximumStepSize = std::numeric_limits::infinity() ; - double InitialStepSize = 1E-8 ; + double InitialStepSize = 1; // Don't make this too small double RelativeTolerance = 1E-15 ; double AbsoluteTolerance = 1E-15 ; @@ -291,10 +266,6 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_VanDerPol ) Eigen::VectorXd Difference = Solution78 - Solution56 ; - std::cout << "Solution = " << Solution78 << std::endl; - std::cout << "Difference between RK56 and RK78: " << std::endl; - std::cout << Difference << std::endl; - BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-13 ) ; BOOST_CHECK_SMALL( std::fabs(Difference(1)) , 1E-13 ) ; } From 2c92a354aa44b1612fc5ae5f78611a392c3d218c Mon Sep 17 00:00:00 2001 From: Rene107 Date: Wed, 23 Mar 2016 08:32:32 +0100 Subject: [PATCH 4/8] Add Logarithmic 2D ODE with analytical solution --- .../numericalIntegratorTestFunctions.h | 66 +++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h b/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h index 9632a6565d..60d42b0471 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h @@ -37,6 +37,7 @@ #define TUDAT_NUMERICAL_INTEGRATOR_TEST_FUNCTIONS_H #include +#include #include @@ -108,6 +109,71 @@ static inline Eigen::VectorXd computeVanDerPolStateDerivative( const double time return vanDerPolStateDerivative; } +//! Computes state derivative of a Logarithmic 2-dimensional ODE +/*! + * Computes the state derivative of a 2-dimensional ODE that is autonomous and includes log functions. + * This function is used by Fehlberg for testing the RK5(6) integrator and has an analytical solution. + * + * \f{eqnarray*}{ + * \frac{ dx }{ dt }( 0 ) &=& -2 * t * x( 0 ) * log( x(1) ) \\ + * \frac{ dx }{ dt }( 0 ) &=& 2 * t * x( 1 ) * log( x(0) ) \\ + * \f} + * + * Reference: Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control + * page 30 + * + * \param time Time at which the state derivative needs to be evaluated. + * \param state State at which the state derivative needs to be evaluated. + * \return Computed state derivative. + */ +static inline Eigen::VectorXd computeFehlbergLogirithmicTestODEStateDerivative( const double time, + const Eigen::VectorXd& state){ + // Declare state derivative vector + Eigen::VectorXd stateDerivative( state.rows() ); + + // Compute state derivative + stateDerivative( 0 ) = -2.0 * time * state(0) * std::log( state(1) ); + stateDerivative( 1 ) = 2.0 * time * state(1) * std::log( state(0) ); + + // Return computed state derivative. + return stateDerivative; +} + +//! Computes the analytical solution of the state for the Logarithmic test ODE of Fehlberg +/*! + * Analytical solution: + * \f{eqnarray*}{ + * x(0) = exp( cos( t^{2} ) ) \\ + * x(1) = exp( sin( t^{2} ) ) \\ + * \f} + * + * with initial conditions: + * \f{eqnarray*}{ + * x_{0}(0) = exp(1) \\ + * x_{0}(1) = 1 + * \f} + * + * Reference: Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control + * page 30 + * + * \param time Time at which the state derivative needs to be evaluated. + * \param initialState State at which the state derivative needs to be evaluated. + * \return Computed state. + */ + +static inline Eigen::VectorXd computeAnalyticalStateFehlbergODE( const double time, + const Eigen::VectorXd& initialState){ + // Declare state derivative vector + Eigen::VectorXd state( initialState.rows() ); + + // Compute state derivative + state(0) = exp( cos( pow(time,2.0) ) ) ; // x(0) = exp( cos( t^{2} ) ) + state(1) = exp( sin( pow(time,2.0) ) ) ; // x(1) = exp( sin( t^{2} ) ) + + // Return computed state derivative. + return state; +} + //! Compute non-autonomous model state derivative. /*! * Computes the state derivative for a non-autonomous ordinary differential equation. The state From 47070a0040fad4023c10cb28dded8e5ce381c47a Mon Sep 17 00:00:00 2001 From: Rene107 Date: Wed, 23 Mar 2016 08:54:21 +0100 Subject: [PATCH 5/8] Changed some comments --- .../UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp index 134b0c00eb..14396e7396 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp @@ -90,7 +90,7 @@ using numerical_integrator_test_functions::computeVanDerPolStateDerivative ; using numerical_integrator_test_functions::computeFehlbergLogirithmicTestODEStateDerivative ; using numerical_integrator_test_functions::computeAnalyticalStateFehlbergODE ; -//! Compare with analytical solution +//! Compare with analytical solution of Fehlberg BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) { tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = @@ -132,8 +132,8 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) Eigen::VectorXd FehlbergError(2) ; FehlbergError << 0.1072E-12 , -0.2190E-12 ; - // sign check -> same sign -> positive - // Not always same sign -> initial step size = 1 or 1E-2, failure computedError(1)/FehlbergError(1) + // sign check + // Not always same sign -> initial step size = 1 or 1E-2, failure: computedError(1) , FehlbergError(1) not same sign BOOST_CHECK_GE(computedError(0)/FehlbergError(0) ,0.0); BOOST_CHECK_GE(computedError(1)/FehlbergError(1) ,0.0); From d51f17abcc8b6f90b501150a0430a67483d9a549 Mon Sep 17 00:00:00 2001 From: Rene107 Date: Wed, 23 Mar 2016 08:56:04 +0100 Subject: [PATCH 6/8] Added a test with a 2D ODE with analytical solution --- ...unitTestRungeKuttaFehlberg78Integrator.cpp | 39 +++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp index 1c57ba3f17..42000ac88b 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp @@ -84,6 +84,45 @@ using numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd; using numerical_integrators::RungeKuttaCoefficients; using numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative; +using numerical_integrator_test_functions::computeFehlbergLogirithmicTestODEStateDerivative ; +using numerical_integrator_test_functions::computeAnalyticalStateFehlbergODE ; + +//! Test Runge-Kutta-Fehlberg 78 integrator using benchmark ODE of Fehlberg(1968) +BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg78_Integrator_Fehlberg_Benchmark ) +{ + tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients::get( + tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + + // Integrator settings + double minimumStepSize = std::numeric_limits::epsilon() ; + double maximumStepSize = std::numeric_limits::infinity() ; + double initialStepSize = 1E-6; // Error: 0.0521 for initialStepSize = 1 ? + double relativeTolerance = 1E-16; + double absoluteTolerance = 1E-16; + + // Initial conditions + double initialTime = 0.0; + double finalTime = 5.0; + Eigen::Vector2d initialState(exp(1.0), 1.0); + + // Setup integrator + tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( + coeff78, computeFehlbergLogirithmicTestODEStateDerivative , + initialTime , initialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); + + + // Obtain numerical solution + Eigen::Vector2d numSol = integrator78.integrateTo(finalTime, initialStepSize); + + // Analytical solution + Eigen::Vector2d anaSol = computeAnalyticalStateFehlbergODE(finalTime, initialState) ; + + Eigen::Vector2d computedError = numSol - anaSol; + BOOST_CHECK_SMALL( std::fabs(computedError(0)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs(computedError(1)) , 1E-13 ) ; +} //! Test Runge-Kutta-Fehlberg 78 integrator using benchmark data from (The MathWorks, 2012). BOOST_AUTO_TEST_CASE( testRungeKuttaFehlberg78IntegratorUsingMatlabData ) From 3b606a31596c994a353c5c8f6c6dbcdaa2a48fea Mon Sep 17 00:00:00 2001 From: Reneh107 Date: Wed, 4 May 2016 12:12:39 +0200 Subject: [PATCH 7/8] Minor style changes --- ...unitTestRungeKuttaFehlberg56Integrator.cpp | 124 +++++++++--------- 1 file changed, 62 insertions(+), 62 deletions(-) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp index 14396e7396..e0a1dc3338 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp @@ -117,29 +117,29 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) // Obtain numerical solution - Eigen::Vector2d numSol = integrator56.integrateTo(finalTime, initialStepSize); + Eigen::Vector2d numericalSolution = integrator56.integrateTo(finalTime, initialStepSize); // Analytical solution // (page 30, Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control) - Eigen::Vector2d anaSol = computeAnalyticalStateFehlbergODE(finalTime, initialState) ; + Eigen::Vector2d analyticalSolution = computeAnalyticalStateFehlbergODE(finalTime, initialState) ; - Eigen::Vector2d computedError = numSol - anaSol; + Eigen::Vector2d computedError = numericalSolution - analyticalSolution; BOOST_CHECK_SMALL( std::fabs(computedError(0)) , 1E-12 ) ; BOOST_CHECK_SMALL( std::fabs(computedError(1)) , 1E-12 ) ; // Error calculated by -> Fehlberg, E. (1968) page 30 // Initial stepsize unknown.. - Eigen::VectorXd FehlbergError(2) ; - FehlbergError << 0.1072E-12 , -0.2190E-12 ; + Eigen::VectorXd fehlbergError(2) ; + fehlbergError << 0.1072E-12 , -0.2190E-12 ; // sign check - // Not always same sign -> initial step size = 1 or 1E-2, failure: computedError(1) , FehlbergError(1) not same sign - BOOST_CHECK_GE(computedError(0)/FehlbergError(0) ,0.0); - BOOST_CHECK_GE(computedError(1)/FehlbergError(1) ,0.0); + // Not always same sign -> initial step size = 1 or 1E-2, failure: computedError(1) , fehlbergError(1) not same sign + BOOST_CHECK_GE(computedError(0)/fehlbergError(0) ,0.0); + BOOST_CHECK_GE(computedError(1)/fehlbergError(1) ,0.0); // Check error is similar - BOOST_CHECK_SMALL(std::fabs(computedError(0)/FehlbergError(0)), 3.0); - BOOST_CHECK_SMALL(std::fabs(computedError(1)/FehlbergError(1)), 1.0); + BOOST_CHECK_SMALL(std::fabs(computedError(0)/fehlbergError(0)), 3.0); + BOOST_CHECK_SMALL(std::fabs(computedError(1)/fehlbergError(1)), 1.0); } @@ -147,127 +147,127 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) { // Setup integrator - tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = tudat::numerical_integrators::RungeKuttaCoefficients::get( tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; - tudat::numerical_integrators::RungeKuttaCoefficients Coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = tudat::numerical_integrators::RungeKuttaCoefficients::get( tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; // Integrator settings - double MinimumStepSize = std::numeric_limits::epsilon() ; - double MaximumStepSize = std::numeric_limits::infinity() ; - double InitialStepSize = 1E-4; // Don't make this too small - double RelativeTolerance = 1E-15 ; - double AbsoluteTolerance = 1E-15 ; + double minimumStepSize = std::numeric_limits::epsilon() ; + double maximumStepSize = std::numeric_limits::infinity() ; + double initialStepSize = 1E-4; // Don't make this too small + double relativeTolerance = 1E-15 ; + double absoluteTolerance = 1E-15 ; // Initial conditions - double InitialTime = 0.5 ; + double initialTime = 0.5 ; Eigen::VectorXd InitialState(1) ; InitialState << 0.5 ; // 1 large error // Setup integrator tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - Coeff56, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, - MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + coeff56, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - Coeff78, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, - MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + coeff78, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); - double EndTime = 1.5 ; - Eigen::VectorXd Solution56 = integrator56.integrateTo(EndTime,InitialStepSize) ; - Eigen::VectorXd Solution78 = integrator78.integrateTo(EndTime,InitialStepSize) ; + double endTime = 1.5 ; + Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize) ; + Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize) ; - Eigen::VectorXd Difference = Solution78 - Solution56 ; + Eigen::VectorXd difference = solution78 - solution56 ; - BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs(difference(0)) , 1E-13 ) ; } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_v2 ) { // Setup integrator - tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = tudat::numerical_integrators::RungeKuttaCoefficients::get( tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; - tudat::numerical_integrators::RungeKuttaCoefficients Coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = tudat::numerical_integrators::RungeKuttaCoefficients::get( tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; // Integrator settings - double MinimumStepSize = std::numeric_limits::epsilon() ; - double MaximumStepSize = std::numeric_limits::infinity() ; - double InitialStepSize = 1; // Don't make this too small - double RelativeTolerance = 1E-10 ; - double AbsoluteTolerance = 1E-10 ; + double minimumStepSize = std::numeric_limits::epsilon() ; + double maximumStepSize = std::numeric_limits::infinity() ; + double initialStepSize = 1.0 ; // Don't make this too small + double relativeTolerance = 1E-10 ; + double absoluteTolerance = 1E-10 ; // Initial conditions - double InitialTime = 0.2 ; + double initialTime = 0.2 ; Eigen::VectorXd InitialState(1) ; InitialState << -1.0 ; // Setup integrator tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - Coeff56, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, - MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + coeff56, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - Coeff78, computeNonAutonomousModelStateDerivative, InitialTime , InitialState, MinimumStepSize, - MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + coeff78, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); - double EndTime = 2.0 ; - Eigen::VectorXd Solution56 = integrator56.integrateTo(EndTime,InitialStepSize) ; - Eigen::VectorXd Solution78 = integrator78.integrateTo(EndTime,InitialStepSize) ; + double endTime = 2.0 ; + Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize) ; + Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize) ; - Eigen::VectorXd Difference = Solution78 - Solution56 ; + Eigen::VectorXd difference = solution78 - solution56 ; - BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-8 ) ; + BOOST_CHECK_SMALL( std::fabs(difference(0)) , 1E-8 ) ; } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_VanDerPol ) { // Setup integrator - tudat::numerical_integrators::RungeKuttaCoefficients Coeff56 = + tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = tudat::numerical_integrators::RungeKuttaCoefficients::get( tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; - tudat::numerical_integrators::RungeKuttaCoefficients Coeff78 = + tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = tudat::numerical_integrators::RungeKuttaCoefficients::get( tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; // Integrator settings - double MinimumStepSize = std::numeric_limits::epsilon() ; - double MaximumStepSize = std::numeric_limits::infinity() ; - double InitialStepSize = 1; // Don't make this too small - double RelativeTolerance = 1E-15 ; - double AbsoluteTolerance = 1E-15 ; + double minimumStepSize = std::numeric_limits::epsilon() ; + double maximumStepSize = std::numeric_limits::infinity() ; + double initialStepSize = 1; // Don't make this too small + double relativeTolerance = 1E-15 ; + double absoluteTolerance = 1E-15 ; // Initial conditions - double InitialTime = 0.2 ; + double initialTime = 0.2 ; Eigen::VectorXd InitialState(2) ; InitialState << -1.0 , 1.0 ; // Setup integrator tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - Coeff56, computeVanDerPolStateDerivative, InitialTime , InitialState, MinimumStepSize, - MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + coeff56, computeVanDerPolStateDerivative, initialTime , InitialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - Coeff78, computeVanDerPolStateDerivative, InitialTime , InitialState, MinimumStepSize, - MaximumStepSize, RelativeTolerance, AbsoluteTolerance ); + coeff78, computeVanDerPolStateDerivative, initialTime , InitialState, minimumStepSize, + maximumStepSize, relativeTolerance, absoluteTolerance ); - double EndTime = 1.4 ; - Eigen::VectorXd Solution56 = integrator56.integrateTo(EndTime,InitialStepSize) ; - Eigen::VectorXd Solution78 = integrator78.integrateTo(EndTime,InitialStepSize) ; + double endTime = 1.4 ; + Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize) ; + Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize) ; - Eigen::VectorXd Difference = Solution78 - Solution56 ; + Eigen::VectorXd difference = solution78 - solution56 ; - BOOST_CHECK_SMALL( std::fabs(Difference(0)) , 1E-13 ) ; - BOOST_CHECK_SMALL( std::fabs(Difference(1)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs(difference(0)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs(difference(1)) , 1E-13 ) ; } BOOST_AUTO_TEST_SUITE_END( ) From aa3576df637b7318ecfd4bee33c6c3f23d5dd1f4 Mon Sep 17 00:00:00 2001 From: Jacco Geul Date: Fri, 13 May 2016 10:03:47 +0200 Subject: [PATCH 8/8] Codecheck --- .../numericalIntegratorTestFunctions.h | 14 +- ...unitTestRungeKuttaFehlberg56Integrator.cpp | 205 +++++++++--------- ...unitTestRungeKuttaFehlberg78Integrator.cpp | 30 +-- 3 files changed, 124 insertions(+), 125 deletions(-) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h b/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h index 60d42b0471..1f3efa6cce 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/numericalIntegratorTestFunctions.h @@ -37,7 +37,6 @@ #define TUDAT_NUMERICAL_INTEGRATOR_TEST_FUNCTIONS_H #include -#include #include @@ -119,15 +118,16 @@ static inline Eigen::VectorXd computeVanDerPolStateDerivative( const double time * \frac{ dx }{ dt }( 0 ) &=& 2 * t * x( 1 ) * log( x(0) ) \\ * \f} * - * Reference: Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control - * page 30 + * Reference: Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta + * Formulas with Stepsize Control (page 30) * * \param time Time at which the state derivative needs to be evaluated. * \param state State at which the state derivative needs to be evaluated. * \return Computed state derivative. */ -static inline Eigen::VectorXd computeFehlbergLogirithmicTestODEStateDerivative( const double time, - const Eigen::VectorXd& state){ +static inline Eigen::VectorXd + computeFehlbergLogirithmicTestODEStateDerivative( const double time, + const Eigen::VectorXd& state){ // Declare state derivative vector Eigen::VectorXd stateDerivative( state.rows() ); @@ -153,8 +153,8 @@ static inline Eigen::VectorXd computeFehlbergLogirithmicTestODEStateDerivative( * x_{0}(1) = 1 * \f} * - * Reference: Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control - * page 30 + * Reference: Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta + * Formulas with Stepsize Control (page 30) * * \param time Time at which the state derivative needs to be evaluated. * \param initialState State at which the state derivative needs to be evaluated. diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp index e0a1dc3338..32c040899e 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp @@ -26,7 +26,7 @@ * YYMMDD Author Comment * 120203 B. Tong Minh Copied RungeKutta4Stepsize unit test. * 120207 K. Kumar Adapted to use modified benchmark functions in Tudat Core. - * 120213 K. Kumar Modified getCurrentInterval() to getIndependentVariable(); + * 120213 K. Kumar Modified getCurrentInterval( ) to getIndependentVariable( ); * transferred to Boost unit test framework. * 120321 K. Kumar Updated (Burden and Faires, 2011) benchmark function call. * 120323 K. Kumar Rewrote unit tests to use benchmark data from @@ -64,11 +64,10 @@ #define BOOST_TEST_MAIN -#include #include #include -#include +#include #include @@ -86,20 +85,21 @@ namespace unit_tests BOOST_AUTO_TEST_SUITE( test_runge_kutta_fehlberg_56_integrator ) using numerical_integrator_test_functions::computeNonAutonomousModelStateDerivative; -using numerical_integrator_test_functions::computeVanDerPolStateDerivative ; -using numerical_integrator_test_functions::computeFehlbergLogirithmicTestODEStateDerivative ; -using numerical_integrator_test_functions::computeAnalyticalStateFehlbergODE ; +using numerical_integrator_test_functions::computeVanDerPolStateDerivative; +using numerical_integrator_test_functions::computeFehlbergLogirithmicTestODEStateDerivative; +using numerical_integrator_test_functions::computeAnalyticalStateFehlbergODE; + +using namespace numerical_integrators; //! Compare with analytical solution of Fehlberg BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) { - tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + RungeKuttaCoefficients coeff56 = + RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg56); // Integrator settings - double minimumStepSize = std::numeric_limits::epsilon() ; - double maximumStepSize = std::numeric_limits::infinity() ; + double minimumStepSize = std::numeric_limits< double >::epsilon( ); + double maximumStepSize = std::numeric_limits< double >::infinity( ); double initialStepSize = 1E-6; // Don't make this too small double relativeTolerance = 1E-16; double absoluteTolerance = 1E-16; @@ -107,39 +107,42 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) // Initial conditions double initialTime = 0.0; double finalTime = 5.0; - Eigen::Vector2d initialState(exp(1.0), 1.0); + Eigen::Vector2d initialState( exp( 1.0 ), 1.0); // Setup integrator - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - coeff56, computeFehlbergLogirithmicTestODEStateDerivative , - initialTime , initialState, minimumStepSize, + RungeKuttaVariableStepSizeIntegratorXd integrator56( + coeff56, computeFehlbergLogirithmicTestODEStateDerivative, + initialTime, initialState, minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); // Obtain numerical solution - Eigen::Vector2d numericalSolution = integrator56.integrateTo(finalTime, initialStepSize); + Eigen::Vector2d numericalSolution = integrator56.integrateTo( finalTime, initialStepSize ); // Analytical solution - // (page 30, Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control) - Eigen::Vector2d analyticalSolution = computeAnalyticalStateFehlbergODE(finalTime, initialState) ; + // (page 30, Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta + // Formulas with Stepsize Control) + Eigen::Vector2d analyticalSolution = + computeAnalyticalStateFehlbergODE( finalTime, initialState ); Eigen::Vector2d computedError = numericalSolution - analyticalSolution; - BOOST_CHECK_SMALL( std::fabs(computedError(0)) , 1E-12 ) ; - BOOST_CHECK_SMALL( std::fabs(computedError(1)) , 1E-12 ) ; + BOOST_CHECK_SMALL( std::fabs( computedError( 0 ) ), 1E-12 ); + BOOST_CHECK_SMALL( std::fabs( computedError( 1 ) ), 1E-12 ); // Error calculated by -> Fehlberg, E. (1968) page 30 // Initial stepsize unknown.. - Eigen::VectorXd fehlbergError(2) ; - fehlbergError << 0.1072E-12 , -0.2190E-12 ; - - // sign check - // Not always same sign -> initial step size = 1 or 1E-2, failure: computedError(1) , fehlbergError(1) not same sign - BOOST_CHECK_GE(computedError(0)/fehlbergError(0) ,0.0); - BOOST_CHECK_GE(computedError(1)/fehlbergError(1) ,0.0); - - // Check error is similar - BOOST_CHECK_SMALL(std::fabs(computedError(0)/fehlbergError(0)), 3.0); - BOOST_CHECK_SMALL(std::fabs(computedError(1)/fehlbergError(1)), 1.0); + Eigen::VectorXd fehlbergError( 2 ); + fehlbergError << 0.1072E-12, -0.2190E-12; + + // Sign check + // Not always same sign -> initial step size = 1 or 1E-2, failure: computedError( 1 ), + // fehlbergError( 1 ) not same sign + BOOST_CHECK_GE( computedError( 0 ) / fehlbergError( 0 ), 0.0 ); + BOOST_CHECK_GE( computedError( 1 ) / fehlbergError( 1 ), 0.0 ); + + // Check error is similar in magnitude + BOOST_CHECK_SMALL( std::fabs( computedError( 0 ) / fehlbergError( 0 ) ), 3.0); + BOOST_CHECK_SMALL( std::fabs( computedError( 1 ) / fehlbergError( 1 ) ), 1.0); } @@ -147,127 +150,123 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Fehlberg_Benchmark ) BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78 ) { // Setup integrator - tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + RungeKuttaCoefficients coeff56 = + RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg56); - tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + RungeKuttaCoefficients coeff78 = + RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg78); // Integrator settings - double minimumStepSize = std::numeric_limits::epsilon() ; - double maximumStepSize = std::numeric_limits::infinity() ; + double minimumStepSize = std::numeric_limits< double >::epsilon( ); + double maximumStepSize = std::numeric_limits< double >::infinity( ); double initialStepSize = 1E-4; // Don't make this too small - double relativeTolerance = 1E-15 ; - double absoluteTolerance = 1E-15 ; + double relativeTolerance = 1E-15; + double absoluteTolerance = 1E-15; // Initial conditions - double initialTime = 0.5 ; - Eigen::VectorXd InitialState(1) ; - InitialState << 0.5 ; // 1 large error + double initialTime = 0.5; + Eigen::VectorXd InitialState( 1 ); + InitialState << 0.5; // 1 large error // Setup integrator - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - coeff56, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, - maximumStepSize, relativeTolerance, absoluteTolerance ); + RungeKuttaVariableStepSizeIntegratorXd integrator56( + coeff56, computeNonAutonomousModelStateDerivative, initialTime, InitialState, + minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - coeff78, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, - maximumStepSize, relativeTolerance, absoluteTolerance ); + RungeKuttaVariableStepSizeIntegratorXd integrator78( + coeff78, computeNonAutonomousModelStateDerivative, initialTime, InitialState, + minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - double endTime = 1.5 ; - Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize) ; - Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize) ; + double endTime = 1.5; + Eigen::VectorXd solution56 = integrator56.integrateTo( endTime, initialStepSize ); + Eigen::VectorXd solution78 = integrator78.integrateTo( endTime, initialStepSize ); - Eigen::VectorXd difference = solution78 - solution56 ; + Eigen::VectorXd difference = solution78 - solution56; - BOOST_CHECK_SMALL( std::fabs(difference(0)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs( difference( 0 ) ), 1E-13 ); } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_v2 ) { // Setup integrator - tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + RungeKuttaCoefficients coeff56 = + RungeKuttaCoefficients::get( + RungeKuttaCoefficients::rungeKuttaFehlberg56 ); - tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + RungeKuttaCoefficients coeff78 = + RungeKuttaCoefficients::get( + RungeKuttaCoefficients::rungeKuttaFehlberg78 ); // Integrator settings - double minimumStepSize = std::numeric_limits::epsilon() ; - double maximumStepSize = std::numeric_limits::infinity() ; - double initialStepSize = 1.0 ; // Don't make this too small - double relativeTolerance = 1E-10 ; - double absoluteTolerance = 1E-10 ; + double minimumStepSize = std::numeric_limits< double >::epsilon( ); + double maximumStepSize = std::numeric_limits< double >::infinity( ); + double initialStepSize = 1.0; // Don't make this too small + double relativeTolerance = 1E-10; + double absoluteTolerance = 1E-10; // Initial conditions - double initialTime = 0.2 ; - Eigen::VectorXd InitialState(1) ; - InitialState << -1.0 ; + double initialTime = 0.2; + Eigen::VectorXd InitialState( 1 ); + InitialState << -1.0; // Setup integrator - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - coeff56, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, - maximumStepSize, relativeTolerance, absoluteTolerance ); + RungeKuttaVariableStepSizeIntegratorXd integrator56( + coeff56, computeNonAutonomousModelStateDerivative, initialTime, InitialState, + minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - coeff78, computeNonAutonomousModelStateDerivative, initialTime , InitialState, minimumStepSize, - maximumStepSize, relativeTolerance, absoluteTolerance ); + RungeKuttaVariableStepSizeIntegratorXd integrator78( + coeff78, computeNonAutonomousModelStateDerivative, initialTime, InitialState, + minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - double endTime = 2.0 ; - Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize) ; - Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize) ; + double endTime = 2.0; + Eigen::VectorXd solution56 = integrator56.integrateTo( endTime, initialStepSize ); + Eigen::VectorXd solution78 = integrator78.integrateTo( endTime, initialStepSize ); - Eigen::VectorXd difference = solution78 - solution56 ; + Eigen::VectorXd difference = solution78 - solution56; - BOOST_CHECK_SMALL( std::fabs(difference(0)) , 1E-8 ) ; + BOOST_CHECK_SMALL( std::fabs( difference( 0 ) ), 1E-8 ); } //! Test Compare with Runge Kutta 78 BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg56_Integrator_Compare78_VanDerPol ) { // Setup integrator - tudat::numerical_integrators::RungeKuttaCoefficients coeff56 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg56) ; + RungeKuttaCoefficients coeff56 = + RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg56 ); - tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + RungeKuttaCoefficients coeff78 = + RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg78 ); // Integrator settings - double minimumStepSize = std::numeric_limits::epsilon() ; - double maximumStepSize = std::numeric_limits::infinity() ; + double minimumStepSize = std::numeric_limits< double >::epsilon( ); + double maximumStepSize = std::numeric_limits< double >::infinity( ); double initialStepSize = 1; // Don't make this too small - double relativeTolerance = 1E-15 ; - double absoluteTolerance = 1E-15 ; + double relativeTolerance = 1E-15; + double absoluteTolerance = 1E-15; // Initial conditions - double initialTime = 0.2 ; - Eigen::VectorXd InitialState(2) ; - InitialState << -1.0 , 1.0 ; + double initialTime = 0.2; + Eigen::VectorXd InitialState( 2 ); + InitialState << -1.0, 1.0; // Setup integrator - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator56( - coeff56, computeVanDerPolStateDerivative, initialTime , InitialState, minimumStepSize, + RungeKuttaVariableStepSizeIntegratorXd integrator56( + coeff56, computeVanDerPolStateDerivative, initialTime, InitialState, minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - coeff78, computeVanDerPolStateDerivative, initialTime , InitialState, minimumStepSize, + RungeKuttaVariableStepSizeIntegratorXd integrator78( + coeff78, computeVanDerPolStateDerivative, initialTime, InitialState, minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); - double endTime = 1.4 ; - Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize) ; - Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize) ; + double endTime = 1.4; + Eigen::VectorXd solution56 = integrator56.integrateTo(endTime,initialStepSize); + Eigen::VectorXd solution78 = integrator78.integrateTo(endTime,initialStepSize); - Eigen::VectorXd difference = solution78 - solution56 ; + Eigen::VectorXd difference = solution78 - solution56; - BOOST_CHECK_SMALL( std::fabs(difference(0)) , 1E-13 ) ; - BOOST_CHECK_SMALL( std::fabs(difference(1)) , 1E-13 ) ; + BOOST_CHECK_SMALL( std::fabs( difference( 0 ) ), 1E-13 ); + BOOST_CHECK_SMALL( std::fabs( difference( 1 ) ), 1E-13 ); } BOOST_AUTO_TEST_SUITE_END( ) diff --git a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp index 42000ac88b..c8e14b011f 100644 --- a/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp +++ b/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg78Integrator.cpp @@ -87,16 +87,16 @@ using numerical_integrator_test_functions::computeNonAutonomousModelStateDerivat using numerical_integrator_test_functions::computeFehlbergLogirithmicTestODEStateDerivative ; using numerical_integrator_test_functions::computeAnalyticalStateFehlbergODE ; -//! Test Runge-Kutta-Fehlberg 78 integrator using benchmark ODE of Fehlberg(1968) +//! Test Runge-Kutta-Fehlberg 78 integrator using benchmark ODE of Fehlberg (1968) BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg78_Integrator_Fehlberg_Benchmark ) { - tudat::numerical_integrators::RungeKuttaCoefficients coeff78 = - tudat::numerical_integrators::RungeKuttaCoefficients::get( - tudat::numerical_integrators::RungeKuttaCoefficients::rungeKuttaFehlberg78) ; + using namespace numerical_integrators; + RungeKuttaCoefficients coeff78 = + RungeKuttaCoefficients::get( RungeKuttaCoefficients::rungeKuttaFehlberg78 ); // Integrator settings - double minimumStepSize = std::numeric_limits::epsilon() ; - double maximumStepSize = std::numeric_limits::infinity() ; + double minimumStepSize = std::numeric_limits::epsilon( ); + double maximumStepSize = std::numeric_limits::infinity( ); double initialStepSize = 1E-6; // Error: 0.0521 for initialStepSize = 1 ? double relativeTolerance = 1E-16; double absoluteTolerance = 1E-16; @@ -104,24 +104,24 @@ BOOST_AUTO_TEST_CASE( test_RungeKuttaFehlberg78_Integrator_Fehlberg_Benchmark ) // Initial conditions double initialTime = 0.0; double finalTime = 5.0; - Eigen::Vector2d initialState(exp(1.0), 1.0); + Eigen::Vector2d initialState( exp( 1.0 ), 1.0); // Setup integrator - tudat::numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( - coeff78, computeFehlbergLogirithmicTestODEStateDerivative , - initialTime , initialState, minimumStepSize, + numerical_integrators::RungeKuttaVariableStepSizeIntegratorXd integrator78( + coeff78, computeFehlbergLogirithmicTestODEStateDerivative, + initialTime, initialState, minimumStepSize, maximumStepSize, relativeTolerance, absoluteTolerance ); // Obtain numerical solution - Eigen::Vector2d numSol = integrator78.integrateTo(finalTime, initialStepSize); + Eigen::Vector2d numericalSolution = integrator78.integrateTo( finalTime, initialStepSize ); // Analytical solution - Eigen::Vector2d anaSol = computeAnalyticalStateFehlbergODE(finalTime, initialState) ; + Eigen::Vector2d analyticalSolution = computeAnalyticalStateFehlbergODE( finalTime, initialState ); - Eigen::Vector2d computedError = numSol - anaSol; - BOOST_CHECK_SMALL( std::fabs(computedError(0)) , 1E-13 ) ; - BOOST_CHECK_SMALL( std::fabs(computedError(1)) , 1E-13 ) ; + Eigen::Vector2d computedError = numericalSolution - analyticalSolution; + BOOST_CHECK_SMALL( std::fabs(computedError( 0 )), 1E-13 ); + BOOST_CHECK_SMALL( std::fabs(computedError( 1 )), 1E-13 ); } //! Test Runge-Kutta-Fehlberg 78 integrator using benchmark data from (The MathWorks, 2012).