# BigFloat promotion rules and constants in functions

**URL:** <https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573>\
**Category:** New to Julia\
**Created:** [September 5, 2018, 10:02am UTC](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573 "2018-09-05T10:02:05Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [September 5, 2018, 11:43pm UTC](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573/21 "2018-09-05T23:43:24Z")

</div>

Try TaylorIntegration.jl for integrators of arbitrarily high order. Cc @lbenet

---

<div class="post-metadata">

**Author:** ![ckoe-bccms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ckoe-bccms/32/1816_2.png) [@ckoe-bccms](https://discourse.julialang.org/u/ckoe-bccms)\
**Post date:** [September 6, 2018, 6:44am UTC](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573/22 "2018-09-06T06:44:16Z")

</div>

> [@platawiec](#):
>
> case, have you tried the [symplectic integrators?](http://docs.juliadiffeq.org/latest/solvers/dynamical_solve.html)

> [@platawiec](#):
>
> s better to use callbacks for energy conservation

Thank you, Henon-Heiles is trivial, I have that working with and without symplectic integrators. The circular restricted three-body problem however is nasty:

The usual “synodic” coordinate system is not an inertial system, it is rotating to seperate the motion of the heavy masses off. The heavy masses are at rest making the ODE for the test particle autonomous and the problem then has only four degrees of freedom. This however has a consequence:

Energy is _not_ an integral of motion ! Also, angular momentum is not an integral of motion. The conserved quantity is the Jacobi Integral which is basically the sum of energy and angular momentum. There exists a second integral (conserved quantity) but that cannot be expressed analytically as far as I am aware. Further, if you look at the formulation as a Newtonian equation of motion for the test particle (so second order differential equations) you see that the force is velocity dependent. As a consequence common symplectic integrators, e.g Verlets or Forrest-Ruth type, make unfortunately no sense as they conserve the wrong quantity and are not applicable anyway due to the velocity dependence of the force 🙂

I originally was planning this for a lecture demonstration (they did that on 1960s computers, it must be simple, right ?) and in fact it is a nice problem where many things fail. However, without very good preparation and introduction … I will probably pass it and what remains is my personal curiosity.

> [@platawiec](#):
>
> t’s better to use callbacks for energy conservation

I am already using a callback to monitor the Jacobi integral and abort if it deviates too much. Eventually I can play with this. My observation is however that small deviations of the integral (1e-12) already for the unstable orbits completelely change the orbit.

Remark: Of course one could work in cartesian coordinates and an inertial system and integrate alle masses (so 18 degrees of freedom). However, near collision orbits still need regularization as far as I read on this.

---

<div class="post-metadata">

**Author:** ![ckoe-bccms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ckoe-bccms/32/1816_2.png) [@ckoe-bccms](https://discourse.julialang.org/u/ckoe-bccms)\
**Post date:** [September 6, 2018, 7:06am UTC](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573/23 "2018-09-06T07:06:19Z")

</div>

> [@Sukera](#):
>
> Still, I’ve dug around a bit and found [this](https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19670016363.pdf) paper from Zurich,

Yes, Eduard Stiefels work. I know that. I gave up early. Probably everything is in there if one can find and specialise it from the general case 🙂

If you are still interested [this](http://ccar.colorado.edu/imd/2015/documents/CRTBP_Handout.pdf) has a executive level summary. Discussion of regularization is missing, though.

---

<div class="post-metadata">

**Author:** ![ckoe-bccms](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ckoe-bccms/32/1816_2.png) [@ckoe-bccms](https://discourse.julialang.org/u/ckoe-bccms)\
**Post date:** [September 6, 2018, 8:08am UTC](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573/24 "2018-09-06T08:08:49Z")

</div>

> [@dpsanders](#):
>
> Try TaylorIntegration.jl for integrators

That method beyond Euler is new to me as a practical one. The general idea is of course simple enough. From a first look there might be a superficial similarity to the recursion approach. However, if floating point rounding error is the actual problem this will not help.

Addendum:

Reading [this article](http://www.maia.ub.es/~angel/taylor/taylor.pdf)  
and looking at the Van der Pol example I am quite sure that NASAs  
recursion method is a very close relative, see [“Power series solutions  
of the Thiele-Burrau regularized planar restricted three-body  
problem”](https://ntrs.nasa.gov/search.jsp?R=19690008692). The article  
also contains an application to the non-planar RTBP. Sections 5.1.3 and  
5.1.4 discuss floating point arithemtics effects. Quote “Both graphics  
show that the error behaviour seem to be dominated by the noise of the  
floating point arithmetic.”, so that appears to be an issue ?

---

<div class="post-metadata">

**Author:** ![PerezHz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/perezhz/32/1750_2.png) [@PerezHz](https://discourse.julialang.org/u/PerezHz)\
**Post date:** [March 29, 2019, 8:27pm UTC](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573/25 "2019-03-29T20:27:49Z")

</div>

I’d just like to mention that there’s an [example](https://perezhz.github.io/TaylorIntegration.jl/latest/common/) in the [TaylorIntegration.jl](https://github.com/PerezHz/TaylorIntegration.jl) docs where the planar RTBP is integrated. Our intention in that example was to show the interoperability between TaylorIntegration.jl and the DifferentialEquations.jl ecosystem. The example shows that TaylorIntegration.jl is able to preserve the Jacobi integral to \sim 10^{-14} (relative), despite having several close encounters with the secondary. Maybe that example can be of some help.

Of course, when dealing with _really_ close encounters (e.g., if the test particle enters the Hill radius of either primary), then regularization helps a lot to preserve the integrals of motion. Some years ago, I used A. Jorba’s [`Taylor`](http://www.maia.ub.es/~angel/taylor/) package (in C) to integrate the planar RTBP and was able to get nice accuracies (i.e., close to the machine epsilon) when dealing with near-collision orbits with the secondary, with Levi-Civita regularization. What I interpret from the quote you posted from Jorba’s paper, is that Taylor’s method is able to reach, at each step, the roundoff accuracy of 64-bit floating point, which is near the best you can ask for with such arithmetic.

Now, when dealing with BigFloats in Julia, with TaylorIntegration.jl you can adjust the parameters `order` and `abstol` in `taylorinteg`, so that you are able to reach the numerical accuracy of BigFloat. There’s a simple (yet, not very user-friendly) [test](https://github.com/PerezHz/TaylorIntegration.jl/blob/master/test/bigfloats.jl) in TaylorIntegration.jl where we actually check this. There, we compute the [complete elliptic integral of the 1st kind](https://en.wikipedia.org/wiki/Elliptic_integral#Complete_elliptic_integral_of_the_first_kind) K(k) to an accuracy of about \sim 10^{-77} using BigFloats, which we then use to compute the solution to the simple pendulum to the same accuracy.

It would be interesting to implement the Thiele-Burrau regularization for the RTBP with TaylorIntegration.jl and BigFloats and see how it performs for your problem; if there’s still some interest I’d be happy to give it a try!

[Previous page](https://discourse.julialang.org/t/bigfloat-promotion-rules-and-constants-in-functions/14573.md?page=1)