I was playing with N-body simulations of a game called Kerbal Space Program, which itself uses the patched conics approximation. I have read that for long term stability it is best to use symplectic integrators and one of the simplest would be Störmer-Verlet method: $$ \ddot{\vec{x}}_n = F(\vec{x}_n) $$ $$ \vec{x}_{n+1} = 2\vec{x}_n - \vec{x}_{n-1} + F(\vec{x}_n){\Delta t}^2 $$ $$ t_{n+1} = t_n + \Delta t $$ I thought that I might add to this method using a similar approach as how this method can be derived, starting with the Taylor expansions: $$ \vec{x}(t+\Delta t) = \sum^N_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^{N+1}) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = \sum^5_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + \sum^5_{n=0} \frac{(-\Delta t)^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^6) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = 2\vec{x}(t) + \Delta t^2 \frac{\partial^2}{\partial t^2}\vec{x}(t) + \frac{\Delta t^4}{12} \frac{\partial^4}{\partial t^2}\vec{x}(t) + O(\Delta t^6) $$ Since the accelerations can be calculated with $F(\vec{x}(t))$ I only had to find an expression for the forth derivative (jounce or snap). For this I used finite difference: $$ \frac{\partial^4}{\partial t^2}\vec{x}(t) = \frac{\ddot{\vec{x}}(t+\Delta t) - 2\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{\Delta t^2} + O(\Delta t^2) $$ $$ \vec{x}(t+\Delta t) = 2\vec{x}(t) - \vec{x}(t-\Delta t) + \Delta t^2 \frac{\ddot{\vec{x}}(t+\Delta t) + 10\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{12} + O(\Delta t^6) $$ This method would be implicit, since for the next position you would need to know the acceleration at that point, however I thought that higher order precision might give it some advantage over Störmer-Verlet. Similar to Störmer-Verlet I assumed that the global error would only drop two orders, such that this method would be a forth order integrator.

Kepler Orbit

I did a simple benchmark of an eccentric Kepler orbit (eccentricity of 0.75 and 20 orbits) I found that this method seems to be a fifth order method and Verlet, as suspected, a second order method. I derived this from the slope in this log-log plot of the time-step size versus the error (distance between numerical and analytical endpoint):

The actual slopes are $4.98\pm 0.01$ for my method and $1.98\pm 0.01$ for Verlet, where the errors are a 95% confidence bound.

Simple Harmonic Oscillator

I also tested this method with a one dimensional simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$, which has the solution $$ x(t-t_0) = \cos(t). $$ To avoid errors in the second required step I used $x(t_0+h)=\cos(h)$ . Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not. I also looked at the error as a function of time-step size and found that this method seemed to be a forth order method and Verlet a second order method when the end position would be at $x_N=x(t_{end})=0$ (maximum speed), but it seemed to be a eighth order method and Verlet a forth order method when the end position would be at $x_N=x(t_{end})=\pm 1$ (minimal speed), see figure below.

The actual slopes are $8.00\pm 0.01$ and $3.998\pm 0.003$ for my method and $3.992\pm 0.005$ and $1.997\pm 0.002$ for Verlet, where the errors are again a 95% confidence bound.

I think this difference in the order of the error is dependent on the local dynamics, but some how for this ODE this error reduces again when speed drops to zero. Therefore it seems that my initial guess was correct and this method is only a forth order accurate.

So the performance of my method depends on where and which ODE is used. Is this a common among other methods? I have had one course about scientific computing, in which a few integration methods where covert, but none of them where symplectic, so I wonder if this integrator is symplectic, since Verlet is said to be symplectic? And does this method have a known name, since it was fairly easy to derive?

PS: In a similar way an eighth order method can be found, by increasing the Taylor expansion by two orders and by using a forth order approximation of the second derivative of the acceleration and a second order approximation of the forth derivative of the acceleration: $$ \vec{x}_{n+1} = 2\vec{x}_{n} - \vec{x}_{n-1} + \Delta t^2 \frac{3\ddot{\vec{x}}_{n+1} + 48\ddot{\vec{x}}_{n+1/2} + 78\ddot{\vec{x}}_{n} + 48\ddot{\vec{x}}_{n-1/2} + 3\ddot{\vec{x}}_{n-1}}{180} + O(\Delta t^8) $$ However I have not yet figured out a good way to find $\ddot{\vec{x}}_{n+1/2}$.

I was playing with N-body simulations of a game called Kerbal Space Program, which itself uses the patched conics approximation. I have read that for long term stability it is best to use symplectic integrators and one of the simplest would be Störmer-Verlet method: $$ \ddot{\vec{x}}_n = F(\vec{x}_n) $$ $$ \vec{x}_{n+1} = 2\vec{x}_n - \vec{x}_{n-1} + F(\vec{x}_n){\Delta t}^2 $$ $$ t_{n+1} = t_n + \Delta t $$ I thought that I might add to this method using a similar approach as how this method can be derived, starting with the Taylor expansions: $$ \vec{x}(t+\Delta t) = \sum^N_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^{N+1}) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = \sum^5_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + \sum^5_{n=0} \frac{(-\Delta t)^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^6) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = 2\vec{x}(t) + \Delta t^2 \frac{\partial^2}{\partial t^2}\vec{x}(t) + \frac{\Delta t^4}{12} \frac{\partial^4}{\partial t^2}\vec{x}(t) + O(\Delta t^6) $$ Since the accelerations can be calculated with $F(\vec{x}(t))$ I only had to find an expression for the forth derivative (jounce or snap). For this I used finite difference: $$ \frac{\partial^4}{\partial t^4}\vec{x}(t) = \frac{\ddot{\vec{x}}(t+\Delta t) - 2\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{\Delta t^2} + O(\Delta t^2) $$ $$ \vec{x}(t+\Delta t) = 2\vec{x}(t) - \vec{x}(t-\Delta t) + \Delta t^2 \frac{\ddot{\vec{x}}(t+\Delta t) + 10\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{12} + O(\Delta t^6) $$ This method would be implicit, since for the next position you would need to know the acceleration at that point, however I thought that higher order precision might give it some advantage over Störmer-Verlet. Similar to Störmer-Verlet I assumed that the global error would only drop two orders, such that this method would be a forth order integrator.

Kepler Orbit

I did a simple benchmark of an eccentric Kepler orbit (eccentricity of 0.75 and 20 orbits) I found that this method seems to be a fifth order method and Verlet, as suspected, a second order method. I derived this from the slope in this log-log plot of the time-step size versus the error (distance between numerical and analytical endpoint):

The actual slopes are $4.98\pm 0.01$ for my method and $1.98\pm 0.01$ for Verlet, where the errors are a 95% confidence bound.

Simple Harmonic Oscillator

I also tested this method with a one dimensional simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$, which has the solution $$ x(t-t_0) = \cos(t). $$ To avoid errors in the second required step I used $x(t_0+h)=\cos(h)$ . Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not. I also looked at the error as a function of time-step size and found that this method seemed to be a forth order method and Verlet a second order method when the end position would be at $x_N=x(t_{end})=0$ (maximum speed), but it seemed to be a eighth order method and Verlet a forth order method when the end position would be at $x_N=x(t_{end})=\pm 1$ (minimal speed), see figure below.

The actual slopes are $8.00\pm 0.01$ and $3.998\pm 0.003$ for my method and $3.992\pm 0.005$ and $1.997\pm 0.002$ for Verlet, where the errors are again a 95% confidence bound.

I think this difference in the order of the error is dependent on the local dynamics, but some how for this ODE this error reduces again when speed drops to zero. Therefore it seems that my initial guess was correct and this method is only a forth order accurate.

So the performance of my method depends on where and which ODE is used. Is this a common among other methods? I have had one course about scientific computing, in which a few integration methods where covert, but none of them where symplectic, so I wonder if this integrator is symplectic, since Verlet is said to be symplectic? And does this method have a known name, since it was fairly easy to derive?

PS: In a similar way an eighth order method can be found, by increasing the Taylor expansion by two orders and by using a forth order approximation of the second derivative of the acceleration and a second order approximation of the forth derivative of the acceleration: $$ \vec{x}_{n+1} = 2\vec{x}_{n} - \vec{x}_{n-1} + \Delta t^2 \frac{3\ddot{\vec{x}}_{n+1} + 48\ddot{\vec{x}}_{n+1/2} + 78\ddot{\vec{x}}_{n} + 48\ddot{\vec{x}}_{n-1/2} + 3\ddot{\vec{x}}_{n-1}}{180} + O(\Delta t^8) $$ However I have not yet figured out a good way to find $\ddot{\vec{x}}_{n+1/2}$.

I was playing with N-body simulations of a game called Kerbal Space Program, which itself uses the patched conics approximation. I have read that for long term stability it is best to use symplectic integrators and one of the simplest would be Störmer-Verlet method: $$ \ddot{\vec{x}}_n = F(\vec{x}_n) $$ $$ \vec{x}_{n+1} = 2\vec{x}_n - \vec{x}_{n-1} + F(\vec{x}_n){\Delta t}^2 $$ $$ t_{n+1} = t_n + \Delta t $$ I thought that I might add to this method using a similar approach as how this method can be derived, starting with the Taylor expansions: $$ \vec{x}(t+\Delta t) = \sum^N_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^{N+1}) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = \sum^5_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + \sum^5_{n=0} \frac{(-\Delta t)^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^6) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = 2\vec{x}(t) + \Delta t^2 \frac{\partial^2}{\partial t^2}\vec{x}(t) + \frac{\Delta t^4}{12} \frac{\partial^4}{\partial t^2}\vec{x}(t) + O(\Delta t^6) $$ Since the accelerations can be calculated with $F(\vec{x}(t))$ I only had to find an expression for the forth derivative (jounce or snap). For this I used finite difference: $$ \frac{\partial^4}{\partial t^2}\vec{x}(t) = \frac{\ddot{\vec{x}}(t+\Delta t) - 2\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{\Delta t^2} + O(\Delta t^2) $$ $$ \vec{x}(t+\Delta t) = 2\vec{x}(t) - \vec{x}(t-\Delta t) + \Delta t^2 \frac{\ddot{\vec{x}}(t+\Delta t) + 10\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{12} + O(\Delta t^6) $$ This method would be implicit, since for the next position you would need to know the acceleration at that point, however I thought that higher order precision might give it some advantage over Störmer-Verlet. Similar to Störmer-Verlet I assumed that the global error would only drop two orders, such that this method would be a forth order integrator.

Kepler problem

I did a simple benchmark of an eccentric Kepler orbit (eccentricity of 0.75 and 20 orbits) I found that this method seems to be a fifth order method and Verlet, as suspected, a second order method. I derived this from the slope in this log-log plot of the time-step size versus the error (distance between numerical and analytical endpoint):

Simple Harmonic Oscillator

I also tested this method with a one dimensional simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$, which has the solution $$ x(t-t_0) = \cos(t). $$ To avoid errors in the second required step I used $x(t_0+h)=\cos(h)$ . Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not. I also looked at the error as a function of time-step size and found that this method seemed to be a forth order method and Verlet a second order method when the end position would be at $x_N=x(t_{end})=0$ (maximum speed), but it seemed to be a eighth order method and Verlet a forth order method when the end position would be at $x_N=x(t_{end})=\pm 1$ (minimal speed), see figure below.

I think this difference in the order of the error is dependent on the local dynamics, but some how for this ODE this error reduces again when speed drops to zero. Therefore it seems that my initial guess was correct and this method is only a forth order accurate.

So the performance of my method depends on where and which ODE is used. Is this a common among other methods? I have had one course about scientific computing, in which a few integration methods where covert, but none of them where symplectic, so I wonder if this integrator is symplectic, since Verlet is said to be symplectic? And does this method have a known name, since it was fairly easy to derive?

PS: In a similar way an eighth order method can be found, by increasing the Taylor expansion by two orders and by using a forth order approximation of the second derivative of the acceleration and a second order approximation of the forth derivative of the acceleration: $$ \vec{x}_{n+1} = 2\vec{x}_{n} - \vec{x}_{n-1} + \Delta t^2 \frac{3\ddot{\vec{x}}_{n+1} + 48\ddot{\vec{x}}_{n+1/2} + 78\ddot{\vec{x}}_{n} + 48\ddot{\vec{x}}_{n-1/2} + 3\ddot{\vec{x}}_{n-1}}{180} + O(\Delta t^8) $$ However I have not yet figured out a good way to find $\ddot{\vec{x}}_{n+1/2}$.

I was playing with N-body simulations of a game called Kerbal Space Program, which itself uses the patched conics approximation. I have read that for long term stability it is best to use symplectic integrators and one of the simplest would be Störmer-Verlet method: $$ \ddot{\vec{x}}_n = F(\vec{x}_n) $$ $$ \vec{x}_{n+1} = 2\vec{x}_n - \vec{x}_{n-1} + F(\vec{x}_n){\Delta t}^2 $$ $$ t_{n+1} = t_n + \Delta t $$ I thought that I might add to this method using a similar approach as how this method can be derived, starting with the Taylor expansions: $$ \vec{x}(t+\Delta t) = \sum^N_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^{N+1}) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = \sum^5_{n=0} \frac{\Delta t^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + \sum^5_{n=0} \frac{(-\Delta t)^n}{n!}\frac{\partial^n}{\partial t^n}\vec{x}(t) + O(\Delta t^6) $$ $$ \vec{x}(t+\Delta t) + \vec{x}(t-\Delta t) = 2\vec{x}(t) + \Delta t^2 \frac{\partial^2}{\partial t^2}\vec{x}(t) + \frac{\Delta t^4}{12} \frac{\partial^4}{\partial t^2}\vec{x}(t) + O(\Delta t^6) $$ Since the accelerations can be calculated with $F(\vec{x}(t))$ I only had to find an expression for the forth derivative (jounce or snap). For this I used finite difference: $$ \frac{\partial^4}{\partial t^2}\vec{x}(t) = \frac{\ddot{\vec{x}}(t+\Delta t) - 2\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{\Delta t^2} + O(\Delta t^2) $$ $$ \vec{x}(t+\Delta t) = 2\vec{x}(t) - \vec{x}(t-\Delta t) + \Delta t^2 \frac{\ddot{\vec{x}}(t+\Delta t) + 10\ddot{\vec{x}}(t) + \ddot{\vec{x}}(t-\Delta t)}{12} + O(\Delta t^6) $$ This method would be implicit, since for the next position you would need to know the acceleration at that point, however I thought that higher order precision might give it some advantage over Störmer-Verlet. Similar to Störmer-Verlet I assumed that the global error would only drop two orders, such that this method would be a forth order integrator.

Kepler Orbit

I did a simple benchmark of an eccentric Kepler orbit (eccentricity of 0.75 and 20 orbits) I found that this method seems to be a fifth order method and Verlet, as suspected, a second order method. I derived this from the slope in this log-log plot of the time-step size versus the error (distance between numerical and analytical endpoint):

The actual slopes are $4.98\pm 0.01$ for my method and $1.98\pm 0.01$ for Verlet, where the errors are a 95% confidence bound.

Simple Harmonic Oscillator

I also tested this method with a one dimensional simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$, which has the solution $$ x(t-t_0) = \cos(t). $$ To avoid errors in the second required step I used $x(t_0+h)=\cos(h)$ . Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not. I also looked at the error as a function of time-step size and found that this method seemed to be a forth order method and Verlet a second order method when the end position would be at $x_N=x(t_{end})=0$ (maximum speed), but it seemed to be a eighth order method and Verlet a forth order method when the end position would be at $x_N=x(t_{end})=\pm 1$ (minimal speed), see figure below.

The actual slopes are $8.00\pm 0.01$ and $3.998\pm 0.003$ for my method and $3.992\pm 0.005$ and $1.997\pm 0.002$ for Verlet, where the errors are again a 95% confidence bound.

I think this difference in the order of the error is dependent on the local dynamics, but some how for this ODE this error reduces again when speed drops to zero. Therefore it seems that my initial guess was correct and this method is only a forth order accurate.

So the performance of my method depends on where and which ODE is used. Is this a common among other methods? I have had one course about scientific computing, in which a few integration methods where covert, but none of them where symplectic, so I wonder if this integrator is symplectic, since Verlet is said to be symplectic? And does this method have a known name, since it was fairly easy to derive?

PS: In a similar way an eighth order method can be found, by increasing the Taylor expansion by two orders and by using a forth order approximation of the second derivative of the acceleration and a second order approximation of the forth derivative of the acceleration: $$ \vec{x}_{n+1} = 2\vec{x}_{n} - \vec{x}_{n-1} + \Delta t^2 \frac{3\ddot{\vec{x}}_{n+1} + 48\ddot{\vec{x}}_{n+1/2} + 78\ddot{\vec{x}}_{n} + 48\ddot{\vec{x}}_{n-1/2} + 3\ddot{\vec{x}}_{n-1}}{180} + O(\Delta t^8) $$ However I have not yet figured out a good way to find $\ddot{\vec{x}}_{n+1/2}$.

However after a simple benchmark of an eccentric Kepler orbit (after 10 orbits) I found that this method is a sixth order integrator. I derived this from the slope in this log-log plot of the time-step size versus the error (difference in distance between numerical and analytical endpoint): This slope does not hold up for big step sizes, since the (relative) error is bigger than which means it is bigger than the radius of the orbit and therefore might become unstable. For very small step sizes my method gets a lot on noise, probably due to the limited accuracy in the initial positions.

However I wonder why my method is a sixth order integrator, does the implicitly mean that the error does not build up over time and therefore the global error is equal to the local error? Will this integrator also be symplectic? And does this method have a known name, since it was fairly easy to derive?

Edit: I also did a quick one dimensional test with a simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$. Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not.

Kepler problem

I did a simple benchmark of an eccentric Kepler orbit (eccentricity of 0.75 and 20 orbits) I found that this method seems to be a fifth order method and Verlet, as suspected, a second order method. I derived this from the slope in this log-log plot of the time-step size versus the error (distance between numerical and analytical endpoint):

Simple Harmonic Oscillator

I also tested this method with a one dimensional simple harmonic oscillator, $\ddot{x}=-x$, with the initial conditions $x(t_0)=1$, $\dot{x}(t_0)=0$, which has the solution $$ x(t-t_0) = \cos(t). $$ To avoid errors in the second required step I used $x(t_0+h)=\cos(h)$ . Even with a time step of only one-forth the period of the oscillation, does the amplitude, after 5000 periods, only drop to 0.973. Would this be enough to say that this method is symplectic, since I have only read about the properties of symplectic integrator, but do not know how to tell whether a method would be symplectic or not. I also looked at the error as a function of time-step size and found that this method seemed to be a forth order method and Verlet a second order method when the end position would be at $x_N=x(t_{end})=0$ (maximum speed), but it seemed to be a eighth order method and Verlet a forth order method when the end position would be at $x_N=x(t_{end})=\pm 1$ (minimal speed), see figure below.

I think this difference in the order of the error is dependent on the local dynamics, but some how for this ODE this error reduces again when speed drops to zero. Therefore it seems that my initial guess was correct and this method is only a forth order accurate.

So the performance of my method depends on where and which ODE is used. Is this a common among other methods? I have had one course about scientific computing, in which a few integration methods where covert, but none of them where symplectic, so I wonder if this integrator is symplectic, since Verlet is said to be symplectic? And does this method have a known name, since it was fairly easy to derive?