I'm trying to implement Backward Euler for a simple physics simulation. My understanding is that Forward Euler 'overapproximates' and adds energy to the system, and so should spiral to infinity, but Backward Euler 'underapproximates' and removes energy to damp to zero.
I have a 2D force field centered at $[\begin{smallmatrix}0\\0\end{smallmatrix}]$, which applies a force tangent to the radial vector. So for some position $\vec{r}=[\begin{smallmatrix}x\\y\end{smallmatrix}]$, the acceleration will be in the direction $\vec{a}=[\begin{smallmatrix}y\\-x\end{smallmatrix}]$. As a matrix equation, I believe the formulation is $\vec{a}=\big[\begin{smallmatrix} 0 & 1\\ -1 & 0 \end{smallmatrix}\big]\vec{r}$. I'll use $B=\big[\begin{smallmatrix} 0 & 1\\ -1 & 0 \end{smallmatrix}\big]$ to denote the matrix here on out.
For Backward Euler, I have the following formulation: \begin{equation} \vec{r}_{t+\Delta{t}}=\vec{r}_t+\Delta{t}\vec{v}_{t+\Delta{t}}\\ \vec{v}_{t+\Delta{t}}=\vec{v}_t+\Delta{t}\vec{a}_{t+\Delta{t}}\\ \end{equation}
Thus, plugging in $\vec{a}_{t+\Delta{t}}$ and solving for $\vec{r}_{t+\Delta{t}}$,
\begin{equation} \vec{v}_{t+\Delta{t}}=\vec{v}_t+\Delta{t}B\vec{r}_{t+\Delta{t}}\\ \vec{r}_{t+\Delta{t}}=\vec{r}_t+\Delta{t}(\vec{v}_t+\Delta{t}B\vec{r}_{t+\Delta{t}})\\ \vec{r}_{t+\Delta{t}}=\vec{r}_t+\Delta{t}\vec{v}_t+\Delta{t}^2B\vec{r}_{t+\Delta{t}}\\ \vec{r}_{t+\Delta{t}}-\Delta{t}^2B\vec{r}_{t+\Delta{t}}=\vec{r}_t+\Delta{t}\vec{v}_t\\ (I-\Delta{t}^2B)\vec{r}_{t+\Delta{t}}=\vec{r}_t+\Delta{t}\vec{v}_t\\ \vec{r}_{t+\Delta{t}}=(I-\Delta{t}^2B)^{-1}(\vec{r}_t+\Delta{t}\vec{v}_t)\\ \end{equation}
I then update the particle's new position and velocity with $\vec{r}_{t+\Delta{t}}$ and $\vec{v}_{t+\Delta{t}}$.
The problem I'm having is that when I implement this, the system gains energy each timestep instead of losing it. I expect the particle to spiral in towards the center, but instead it spirals out to infinity just like Forward Euler. I've been trying to search for why this happens but can't seem to figure it out.
Is there something wrong with my formulation? And why is the system not being damped/losing energy?
Thanks for reading.
In complexified coordinates $z=x+yi$ you want to solve $$\ddot z=-iz.$$ This has basis solutions $$ \exp\left(\pm\frac{1+i}{\sqrt2}t\right), $$ so already the exact solution will spiral out, it is not astonishing that a numerical solution follows this pattern.
A radial vector field is always proportional to the radius vector, your vector field is normal to the radius vector. It is tangent to circles with center at the origin.
Your vector field is not conservative, so any argument with changing vs. constant energy levels is invalid.
Always use a better integration method, a higher order method or smaller step sizes, if you think that your implementation has an error. It is then easy to see if the difference between the two is completely out-of-line or falls into an expected error growth pattern.