# How to properly implement Backwards Euler on a system of bodies

I have a system of two bodies, $b^i$ and $b^j$, each in position $\vec{p}^i = (x^i, y^i)$ and $\vec{p}^j = (x^j, y^j)$. The two are connected via spring and I'd like to know, given their states at time $t_n$, their states at time $t_{n+1}$ ($\vec{p}^i_{n+1}$ and $\vec{p}^j_{n+1}$).

I have decided to use the backward Euler method since my spring is going to be stiff. However, I am not sure how to implement it.

Currently, I am implementing it the following way:

$\vec{p}^i_{{n+1}} = \vec{p}^i_{n} + hf(t_n+h, \vec{p}^i_{n+1}, \vec{p}^j_{n})$

$\vec{p}^j_{{n+1}} = \vec{p}^j_{n} + hf(t_n+h, \vec{p}^j_{n+1}, \vec{p}^i_{n})$

As you can see, I am using the other body's state at each backward Euler step. However, my system won't stabilize sometimes.

Am I supposed to, instead, build a matrix $P$

$P = \begin{bmatrix} x^i & y^i \\ x^j & y^j \end{bmatrix}$

an then use it to build a backward Euler step of the form

$P_{n+1} = P_n + hf(t_{n}+h, P_{n+1})$?

If so, how would I implement Newton's method for this matrix? to use Newton's method for a vector of size 2, I build a 2x2 Jacobian matrix, but how do I do it for a matrix of 2x2, do I make some sort of 4x4 Jacobian matrix?

EDIT:

To further elaborate on my function f

$f(t_n, \vec{p}^i_{n}, \vec{p}^j_{n}) = K \times (\vec{p}^j_n - \vec{p}^i_n)$

where $K$ is an attraction coefficient which is user adjustable.

So in this case, f requires both the position of $b^i$ and $b^j$ to evaluate a displacement for each of the bodies

• Backward Euler is easy to implement, if f(x,t) is a simple function. I found one simple example while google it, here is the link. math.la.asu.edu/~dajones/class/275/ch2.pdf – AGN Sep 6 '16 at 4:16
• Did u check whether there is any numerical stability issues in that system for explicit methods? I believe you didn't get confused with spring stiffness and stiff-ODEs – AGN Sep 6 '16 at 4:20
• Yes, there are numerical instability issues with the system if I use explicit Euler, even explicit 4th order Runge Kutta method loses stability if my time-steps are too big. Yes, when I tagged the question with "stiffness" I though this talks about systems with stiff springs. The example you posted is very simple since it contains the computation for only one variable and only one particle, I don't even need a Jacobian matrix for that, I'm interested in solving an expression with multiple particles and multiple coordinates for each, so the example doesn't easily translate – WhiteTiger Sep 6 '16 at 4:49
• For system of equation, we can use the same procedure as long as the system is not coupled, (coupled system need little bit extra maths because can't directly apply linear algebra ). In that example $p$ is scalar, here its a matrix. U can subs. $p^{n+1}$ in f(p,t) and, u can convert that to explicit function like $p^{n+1} =A*p^n+B$ using symbolic tools. Using matlab, mathematica or sympy u can find And B in that expression, then coding is straight forward like explicit schemes – AGN Sep 6 '16 at 5:06
• So you mean I can compute the $P_{n+1}$ inside f(t, p) using Forward Euler and then solve the Backward Euler? I suspect my system IS coupled since, as I mentioned, the forces acting on body $b^i$ depend both on the position of $b^i$ and $b^j$. Also, I need to code this in Java, I'm building a simple physics system for an application so this has to get solved while I'm computing – WhiteTiger Sep 6 '16 at 5:37

I assume

$p=[p^i;p^j]$

$p_i=[x^i, y^i]$

$p_j=[x^j, y^j]$

Formula for backward Euler is:

$p^{n+1}=p^n+h*f(p^{n+1})$

here $f_i=K*(p_i-p_j)$

$f_j=K*(p_j-p_i)$

Substituting this in backward Euler's formula

$$[p_i^{n+1};p_j^{n+1}] =[p_i;p_j]^n+h*K*[p_i^{n+1}-p^n_j;p_j^{n+1}-p^n_i]$$

Implicated only some part of second term in RHS to enhance stability!

Using matrix assoative property, we can write $$[p_i^{n+1};p_j^{n+1}] =[p_i;p_j]^n+h*K*[p_i^{n+1};p_j^{n+1}]-h*K*[p_j^{n};p_i^{n}]$$

$$[p_i^{n+1};p_j^{n+1}](1-Kh)=[p_i;p_j]^n-Kh*[p_j^{n};p_i^{n}]$$

This form looks like explicit schemes and easy to implement. You can directly use this scheme using matrix libraries available in programming languages. For obtaining $p^{n+1} =A*g(p^n)$ foam, you can use symbolic languages like sympy, matlab, mathematica, maple etc, Other wise you can use linear or non-linear root-finding method to solve $p^{n+1}=p^n+h*f(p^{n+1})$ equation iteratively. This is why implicit methods are computationally expensive than explicit method. If we do some little algebra like what we did, you can save computational cost as well as easy to parallelize the code.

Limitations

1. It is a first order scheme may deviate from actual solution.
2. Need fine mesh to get accurate solution.
3. You may consider higher order time integration.

Note: Actual spring mass system we balance inertial force by spring force and inertial force is second order term, for that you can split the second order equation into two first order equation and solve it. I'm not sure whether it is completely implicit or semi-implicit scheme. Semi-implicit scheme is conditionally stable for some problems.

• Okay, I'm confused about this notation a bit. So you wrote $f(P_{n+1}) = K * P_{n+1}$, but that isn't what I wrote as my implementation of $f$. My implementation relies on both $p^i$ (position of body $b^i$) and $p^j$ (position of body $b^j$). So if when you wrote $P_{n+1}$ you meant THE MATRIX $P = [[x^i, y^i], [x^j, y^j]]$ then I'd have $f(P_{n+1}) = K * P_{n+1}[1] - P_{n+1}[0]$ (the position of $b^j$ minus the position of $b^i$). What I'm basically asking is "How do I implement implicit Euler on a system of bodies whose movements are governed by higher order expressions?" – WhiteTiger Sep 6 '16 at 20:02
• @WhiteTiger please check it! – AGN Sep 7 '16 at 3:59