I've read an article (Long-term integrations and stability of planetary orbits in our Solar system) in which the authors solved the problem of the absence of an analytical solution for the solar system planetary dinamics by making a very accurate solution. They called it a "pseudo-exact solution". They then used this solution to estimate the error in planetary longitudes of their main solution, which was less precise. For instance the time step of the "exact" solution was 0.125 days, while the time step of the main solution was 8 days, meaning that the latter was 64 times bigger than the one for "exact" solution.

What I don't understand is how they could compare the two solutions since the time steps was so different. I mean, maybe they did something like (u - y(1:64:N))/(2^p - 1)? Where u is the main solution, y the precise solution and N is the number of iterations. But in that case how they could know that the terms of y to compare with u was precisely every 64 steps? I mean, the second element of u for example could be compared with a term between the index 2 and 64 of y actually. How can be chosen a specific value between 2 and 64 without checking the behaviour of the exact solution in every 64-long interval?

  • 1
    $\begingroup$ if you set u(1) = y(1) and you set your time steps such that h_u = h_y/64, then you expect u(2) ~ y(65). It makes no sense to compare the solutions at different times $\endgroup$ – Yann Feb 19 at 15:04
  • $\begingroup$ @Yann That was just an example, and it was based on this answer: math.stackexchange.com/questions/3058387/… (here instead of $h_y / 64$, $h_y/2$ is used). I know it could not have sense in the case of my question, and that's also why I was wondering how they could have compared that solutions. $\endgroup$ – Zebx Feb 19 at 15:43
  • $\begingroup$ Perhaps the Richardson extrapolation techniques? $\endgroup$ – Maxim Umansky Feb 19 at 16:19
  • $\begingroup$ @MaximUmansky That's what I supposed, even if they used a time step $64$ times bigger than the "exact" solution. Maybe it was not very clear from my example, I edited the question. $\endgroup$ – Zebx Feb 19 at 16:31
  • $\begingroup$ You typically would interpolate the computed solution to get a solution at any desired time. With adaptive step sizes, there is no practical way to ensure that the integrator computes the solution to a desired point in time exactly. $\endgroup$ – Brian Borchers Feb 19 at 21:24

Browsing the paper I can confirm that all integrations of the trajectories is done with a fixed time step.

The authors compute approximations $A_h(t)$ for all $t \in \Sigma_h$ where $$\Sigma_h = \{ nh \: : \: 0 \leq nh \leq T\}.$$ They use at least two different values of $h$, namely $h_1=0.125$ days and $h_2 = 8$ days. It is important to recognize that $$ \Sigma_{h_2} \subset \Sigma_{h_1}.$$ This implies that for all $t \in \Sigma_{h_2}$ we have two different approximations of all target values, i.e., the planetary positions. It is common to compare these values, but it also somewhat dangerous. The problem is that they do not know that the computed values for $h=h_1$ are reliable. It is entirely possible that $h=h_1$ is so small that rounding errors rather than truncation error are dominating for this calculation. I am not saying that this is the case, I am simply observing that this is a possibility.

A safer approach is to compute multiple approximations using time step $h_j = 2^{-j}h_0$ for some suitable value of $h_0$ and $j=0,1,2, \dots.$ If the true value $S(t)$ and the approximation $A_h(t)$ satisfies an asymptotic error expansion of the form $$ S(t) - A_h(t) = \alpha(t) h^p + \beta(t) h^q + O(h^s), \quad 0<p<q<s$$ where $\alpha$ and $\beta$ are functions that are independent of $h$, then Richardson's fraction given by $$ F_h(t) = \frac{A_{2h}(t) - A_{4h}(t)}{A_h(t) - A_{2h}(t)}, $$ will behave in a very predictable manner. Specifically, we have $$ F_h \rightarrow 2^p, \quad h \rightarrow 0, \quad h>0$$ We also have $$ F_h - 2^p = O(h^{q-p}), \quad h \rightarrow 0, \quad h > 0$$ and $$ S(t) - A_h(t) = \frac{A_h(t) - A_{2h}(t)}{2^p-1} + O(h^q), \quad h \rightarrow 0, \quad h > 0.$$ This is true in exact arithmetic, but is also very nearly true in floating point arithmetic up right up to the point where $h$ is so small that rounding errors start to become comparable to the truncation error, i.e., the value $S(t) - A_h(t)$ that we would dearly like to estimate. By tracking the development of Richardson's fraction we can manually determine if we are in the so-called asymptotic range where rounding errors are irrelevant and Richardson's error estimate $$ S(t) - A_h(t) \approx \frac{A_h(t) - A_{2h}(t)}{2^p-1}$$ is a good approximation. In the case of the paper, the authors uses a method that is second order accurate. This corresponds to $p=2$. For such a smooth problem, I would be very surprised if the error does not obey an asymptotic error expansion.

  • $\begingroup$ Yes, regarding the first part of your answer, that was another thing that wasn't clear to me. For instance how could they know their 0.125 days solution was actually good enough. The rest of your answer explain clearly the approach I (poorly) supposed they could have at least used in some way. Now things seems really much clearer to me, thank you very much. A question: how much precise the $2^p$ factor should be? I mean how much could be the deviation of $F_h$ from $2^p$ for Richardson error to still be a good approximation? $\endgroup$ – Zebx Feb 24 at 19:43
  • 1
    $\begingroup$ @Zebx Thank you for your kind words. Your new question is good and it can be answered in terms of the asymptotic error expansion. This is just a sketch. We need the exact formulas for $F_h$ and Richardson's error estimates in terms of $\alpha$, and $\beta$. They will reveal that the behavior of $F_h $ is controlled by $\frac{\beta}{\alpha}h^{q-p}$ which is precisely the term needed to judge the quality of Richardson's error estimate. It is a little bit tedious to write out, but it is essentially just a matter of tracking power of $h$. $\endgroup$ – Carl Christian Feb 24 at 21:04
  • $\begingroup$ Ok, so one can have an idea at least looking at $F_h − 2^p$, considering that $\alpha$ and $\beta$ are not easy to evaluate I suppose. I like this procedure, I'm working on the final project of an university course but we haven't been taught how to estimate numerical error without an analytical solution (yes, it's absurd). I get back to work, let's see what I can get. Thank you again. $\endgroup$ – Zebx Feb 25 at 17:10
  • $\begingroup$ Hi, sorry if I comment here after so much time, could you provide me with some reference about the richardson fraction and the error estimation as you explained in your answer? $\endgroup$ – Zebx Apr 9 at 17:45
  • $\begingroup$ @Zebx: No worries. I am writing a textbook for my own students. If you send me an email, then you can have the current draft of the relevant chapters. $\endgroup$ – Carl Christian Apr 10 at 9:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.