Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Runge Kutta 56 integrator to coefficient list + Unit test #22

Merged
merged 10 commits into from
May 13, 2016

Conversation

Reneh107
Copy link
Contributor

A unit test is written that performs three tests comparing the results of the RK78 and the RK56 integrator.

The RK56 integrator is added to the lists such that it can be included.

Code check required...

Note: Solution was sensitive to integrator settings. Possible due to unstable behavior of ODEs.

Add test case from Fehlberg1968.
@magnific0
Copy link
Member

Hi @Reneh107, thanks a lot for the pull request. The RKF56 has been "disabled" because it lacks a proper unit test. With your work we come a long way. I do have a couple of remarks:

  1. The test does not verify the implementation of RKF56, it proves that a good enough propagator has been implemented. This was the original reason for not enabling 56. I came by a report of Fehlberg hoping for some numerical values that would prove this. See my pull request on your branch and (page 30, Fehlberg, E. (1968). Classical Fifth-, Sixth-, Seventh- and Eigth-Order Runge-Kutta Formalas with Stepsize Control). The prove will rely on reproducing the same error wrt the exact solution under certain settings. However, my results so far have been far off. It would be really great if we can get such a test in there.

  2. The first of the unit tests fails for me, I get the following error message

    Test 1 : Comparison RK56 and RK78
    Solution = 2.52349
    Difference between RK56 and RK78: 
    -0.328654
    [...]/tudat/Tudat/Mathematics/NumericalIntegrators/UnitTests/unitTestRungeKuttaFehlberg56Integrator.cpp(186): error: in "test_runge_kutta_fehlberg_56_integrator/test_RungeKuttaFehlberg56_Integrator_Compare78": absolute value of std::fabs(Difference(0)){0.32865413740151306} exceeds 1e-14
    
  3. Personally I think this is very minor, but it matters nonetheless. The code you have written is not entirely conform the style guide. Most importantly the variable naming convention is lowerUpperCase. Also you seem to have few "[space][semi colon]" here and there. Also consider less verbose unit tests, they're great for debugging, but the convention is to not have std::cout's in the code. Same thing goes for portions of commented code.

Anyway, thanks again. I hope you find these useful. Looking forward to your reaction.

@Reneh107
Copy link
Contributor Author

I processed the comments and updated the unit test.
I added the logarithmic test function to the test functions header.

It was observed that the accuracy of the RK56 numerical integrator depends on the initial stepsize. For some initial stepsizes, a large error was observed, In the order of 0.01. This could mean that there is something wrong with the Error control. However, it could also be an inherent problem of the RK56 integrator if these initial stepsizes correspond to the unstable region of the RK56 integrator.

The same test was also added to the RK78 integrator and it was observed that the accuracy also depends on the initial stepsize.

RK5(6):
Initial stepsize: 1E-10 , Error: 0.0054
Initial stepsize: 1E-6 , Error: <1E-13
Initial stepsize: 1 , Error: <1E-13

RK7(8):
Initial stepsize: 1E-6 , Error: <1E-13
Initial stepsize: 1 , Error: 0.0521

Conclusion:

  1. Integrator works
  2. Error/stepsize control is not optimal OR instability of the RK integrators was observed. (This applies to all RK variable stepsize integrators)

@magnific0 magnific0 assigned magnific0 and unassigned magnific0 May 4, 2016
@Reneh107
Copy link
Contributor Author

Reneh107 commented May 4, 2016

I updated the unit test such that it agrees with the tudat coding style.

I propose that (after code check) we merge this pull request into tudat and open an issue to improve the unit test to verify that the integrator is indeed a RK56 integrator.

The current unit test only proves that the integrator works correctly. Tests have shown that the obtained stepsize with this integrator is in between the stepsize that the 45 and the 78 generate. Therefore, one can assume that this integrator is indeed a correctly implemented 56 integrator.

@magnific0
Copy link
Member

magnific0 commented May 13, 2016

@Reneh107

Thanks for making an effort. This code will be very useful for others and unit testing wise this is as complete as we'll get it now. The idea of opening an issue about this is indeed much more practical.

I went through your code and made a few changes:

  • Styling mostly to do with spacing. Don't precede "," and ";" with spaces. Put space everywhere else :)
  • Watch out for the 100 line width
  • Namespace: you are already in namespace tudat and tudat::unit_tests, so it's not necessary to precede everything with tudat::. Similarly a simple using namespace numerical integrators can make the code much more readable.
  • Last but not least try to include modern implementations of C libraries, so <cmath> instead of <math.h>, rule of thumb: if you include a standard (Tudat and Boost are a different story) library with a .h, there is likely a different package you should use.

I have just made a pull request on your branch. Please accept it such that the changes are reflected here. Then I will merge the code into Tudat!

@Reneh107
Copy link
Contributor Author

Reneh107 commented May 13, 2016

@magnific0

Thank you for the feedback! I will use this in my upcoming pull requests.

I merged your changes into my branch, so you can merge it in Tudat. I will open an issue after you've merged the code.

@magnific0 magnific0 merged commit efd07f2 into Tudat:master May 13, 2016
@gviavattene gviavattene mentioned this pull request Sep 21, 2017
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants