I want to solve a two-body problem using RK2. The following is a system of equations to be solved
\begin{equation}\label{two-body}
x^{\prime \prime}=-\frac{x}{r^{3}}, y^{\prime \prime}=-\frac{y}{r^{3}}, r=\sqrt{x^{2}+y^{2}}
\end{equation}
let
\begin{equation}
\begin{aligned}
v_{1}=x \\
v_{2}=y \\
v_{3}=x^{\prime} \\
v_{4}=y^{\prime}
\end{aligned}
\end{equation}
So we get
\begin{equation}
\begin{aligned}
v_{1}^{\prime}=v_{3} \\
v_{2}^{\prime}=v_{4} \\
v_{3}^{\prime}=-\frac{v_{1}}{r^{3}} \\
v_{4}^{\prime}=-\frac{v_{2}}{r^{3}}
\end{aligned}
\end{equation}
Then let $\textbf{v}=\left[v_1,v_2,v_3,v_4\right]=\left[x, y, x^{\prime}, y^{\prime}\right]$ and $f(t, v)=v^{\prime}(t)=\left[x^{\prime}, y^{\prime},-\frac{x}{r^{3}},-\frac{y}{r^{3}}\right], r=\sqrt{x^{2}+y^{2}}$.
Next, I arrange IVP like this
\begin{equation}\label{MSA two body}
\begin{aligned}
\textbf{v}^{\prime}(t)=f(t, \textbf{v})\\
\textbf{v}\left(t_{0}\right)=\textbf{v}_{0}.
\end{aligned}
\end{equation}
Now, I want to solve with RK2
$
\begin{aligned}
\mathbf{v}_{n+1}=& \mathbf{v}_{n}+h\left(b_{1} k_{1}+b_{2} k_{2}\right) \\
k_{1}=& f\left(t_{n}, \mathbf{v}_{n}\right)=\left[x_{n}^{\prime}, y_{n}^{\prime}, x_{n}^{\prime \prime}, y_{n}^{\prime \prime}\right] \\
k_{2}=& f\left(t_{n}+a_{21} h_{1}, \mathbf{v}_{n}+a_{21} h k_{1}\right) \\
=& {\left[x_{n}^{\prime}+a_{21} h\left(\frac{-x_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}\right), y_{n}^{\prime}+a_{21} h\left(\frac{-y_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}\right),\right.} \\
&\left.\frac{-\left(x_{n}+a_{21} h x_{n}^{\prime}\right)}{\left(\left(x_{n}+a_{21} h x_{n}^{\prime}\right)^{2}+\left(y_{n}+a_{21} h y_{n}^{\prime}\right)^{2}\right)^{\frac{3}{2}}}, \frac{-\left(y_{n}+a_{21} h y_{n}^{\prime}\right)}{\left(\left(x_{n}+a_{21} h x_{n}^{\prime}\right)^{2}+\left(y_{n}+a_{21} h y_{n}^{\prime}\right)^{2}\right)^{\frac{3}{2}}}\right] \\
=& {\left[k_{1 x}+a_{21} h k_{1 x^{\prime}}, k_{1 y}+a_{21} h k_{1 y^{\prime}},\right.} \\
&\left.f_{x^{\prime}}\left(x_{n}+a_{21} h k_{1 x}, y+a_{21} h k_{1 y}\right), f_{y^{\prime}}\left(x_{n}+a_{21} h k_{1 x}, y_{n}+a_{21} h k_{1 y}\right)\right] .
\end{aligned}
$
with $k_{1 x}= x_{n}, k_{1 y}= y_{n}, k_{1 x^{\prime}}= \frac{-x_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}, k_{1 y^{\prime}}= \frac{-y_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}$ and $f_{x^{\prime}}\left(x_{n}, y_{n}\right)= \frac{-x_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}},f_{y^{\prime}}\left(x_{n}, y_{n}\right)= \frac{-y_{n}}{\left(x_{n}^{2}+y_{n}^{2}\right)^{\frac{3}{2}}}$.
Now, I will compute analytics and numeric value (with RK2 there are Heun, Ralston,midpoint and new method with unusual RK2 coefficient) from two-body problem then compare the error. Here the program
https://colab.research.google.com/drive/1P0aRbqmFa1b77mGtzsojmZwe3BThzf_D?usp=sharing
And I got plot like this
The question is how can when h decreasing the error gets bigger in all methods after N>600? (N is number of point on interval)
Best Answer
Two conceptual errors in the code
you start every method step from the exact solution. However, to explore the global error of the methods, you need to continue from the numerical values so that the errors can accumulate.
when you compute the $L^2$ norm of the error, then you need to multiply the square sum by $h$, and return the square root of that, $\|f\|_{L^2}\simeq \sqrt{\sum |f(x_k)|^2\Delta x}$.
In the second point the effects of the first point will distort the picture somewhat. The global error is (weakly) proportional to the unit step error. The single-step error that you compute is $h$ times the unit step error. Thus in the square sum you get a factor $h$ too many.
But in the main it is as said in the comments, the global error behaves similar to $\frac{\mu}{h}+h^2$, $\mu$ the general error per floating point operation, accumulated over $n\sim \frac1h$ steps, and the method step error $\sim h^3$ similarly accumulated to a global error proportional to $h^2$. This gives the typical hook shape with the apex around $h=10^{-5}$ for 64bit double precision.