Comparing numerical solutions with very different time grids

Question

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?

Answer 1

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.