Open In Colab

Remark:

This article is not yet readable. We've hit a technical difficulty which need to first solve, namely how to integrate ODEs of the form $x''=f(x,x', x'')$ on python. Typically odeint integrates equations of the form $x''=f(x)$ or $x''=f(x,x')$. But the Weber force involves interaction between all the particles their relative accelerations, so applying Newton's 2nd Law gives us equations of motion of the form $$m_i a_i = m_i v'_i=\sum_j F_{ij}$$ where the force depends on the acceleration of $a_i$ and the particles $a_j$. In the literature this is known as a "Differential Algebraic Equation" (DAE). So we've been looking for python packages which have DAE solvers. It does not appear to be well-known subject.

Introduction: Weber's Force.

Our goal is to study the $N$-body problem on python according to Weber's electrodynamic law. In fact $N=3$ would be good start, and many author's have studied the $2$-body problem, even us in our first post on this blog.

The Weber potential belongs to relational mechanics, being given by a force $F_{12}$ between two charged particles satisfying Newton's action-reaction principle, so $F_{12}=-F_{21}$, and the force is also radial being along the straight line segment joining the centre of the particles. However what is further strange about Weber's potential is that the force $F$ depends on the relative positions, velocity, and accelerations of the particles!

So let us consider a three body system with particles having inertial masses $m_1, m_2, m_3$ and such that they have positions $r_1, r_2, r_3$ relative to the centre of mass $O=O_{CM}$ of the system. We will assume that the calculations are performed in the centre of mass frame.

The total force acting on particle $m_1$ is the sum $F_{12}+F_{13}$, where for example $$F_{1j}=\frac{q_1 q_j}{r^3_{1j}} \hat{\bf{r}}_ \left( 1-\frac^2}{2 c^2} + \frac{r_{1j} r_{1j}''}{c^2} \right) $$ for $j=2,3$.

Now if we directly apply Newton's 2nd Law we get: $$F_{\text{net} 1}=F_{21}+F_{31}=m_1 \frac{dv_1}{dt}.$$ Here we want to be extremely careful in interpreting these terms. The velocity $v_1$ is the velocity of the particle $m_1$ with respect to the centre-of-mass, i.e. we have $v_1=v_{1O}$. The difficulty with the second order differential equation obtained from Newton's 2nd Law is that $r_{12}''$ is an acceleration term that has some component in common with $\frac{dv_1}{dt}$. The term $r_{12}''$ is very different from $r_1$. Indeed we view $r_1$ and $dr_1/dt=v_1$ as vectors (a position and velocity vector). However the radial $r_{12}$ and $r_{12}'$, $r_{12}''$ are scalar-valued functions!

To use python's odeint solver, we either need to rewrite the expression defining Newton's 2nd Law, or we need to directly define an integration method which integrates the equation as is.

Formulas and Computation:

What are the terms that enter into the ODE?

In this setting a state has the form $$\{(x_i, v_i, m_i, q_i)\}_i$$ for a finite collection of particles having position $x_i$, velocity $v_i$, and mass $m_i$ and electric charge $q_i$. In our computations it's convenient to separate the spatial variables $[x_i, v_i]$ from the mass charge $[m_i, q_i]$ for each $i^{th}$ particle in the system.

We suppose our system is isolated, and therefore the net force acting on the $i^{th}$ particle is the sum of the forces $F_{ij}$ over all indexes $j$. The pairwise particle force $F_{ij}$ is evaluated by Weber's force law. The basic notation is:

  • ${\bf{r}}_{ij} := r_i - r_j ~~~......... ~\text{vector}$
  • $r_{ij}:=|r_i - r_j| ~~~......... ~\text{scalar}$
  • $v_{ij}:=v_i - v_j ~~~......... ~\text{vector}$
  • $a_{ij}:=a_i -a_j ~~~......... ~\text{vector}$
  • $r'_{ij}:= {\hat{\bf{r}}}_{ij} \cdot v_{ij} ~~~......... ~\text{scalar}$
  • $r_{ij}~r''_{ij}:={v_{ij} \cdot v_{ij}} - ({\hat{r}}_{ij} \cdot v_{ij})^2 +{\bf{r_{ij}}} \cdot a_{ij}. ~~~......... ~\text{scalar}$

This last identity for $r''_{ij}$ is important for our calculations. If we directly substitute this formula into Weber's force law and Newton's 2nd Law, then the net force on the $i^{th}$ particle yields:

$$m_i a_i= f_i(r_{ij}, v_{ij}) + \sum_{j} \frac{q_i q_j}{r_{ij}^2 c^2} ~ {\hat{\bf{r}}}_{ij}~ ( r_{ij} \cdot a_{ij}), $$

where $f_i$ is defined as:

$$f_i(\{r_{ij}, v_{ij}\}):=\sum_j \frac{q_i q_j}{c^2 r^2_{ij}} {\hat{\bf{r}}}_{ij} \left( c^2-\frac{3}{2}({\hat{\bf{r}}}_{ij} \cdot v_{ij})^2 + |v_{ij}|^2 \right) $$

for $1\leq i , j \leq 3$.

For every $i$, we want to isolate the acceleration terms $$m_i a_i -\sum_{j} \frac{q_i q_j}{c^2r_{ij}^2} ~ {\hat{\bf{r}}}_{ij}~ ( r_{ij} \cdot a_{ij}) $$ as much as possible, and this leads to our solving Implicit ODEs or so-called DAEs. Ideally we would avoid the DAE solvers as much as possible, and would prefer to have all acceleration terms $\{a_i\}_i$ isolated. But it's not clear a priori whether this is possible.

Using the above formulas we rewrite $r'_{ij} = \hat{\bf{r}}_{ij} \cdot v_{ij}$, and therefore everything in $f_i$ can be computed directly from $ \{ {\hat{\bf{r}}}_{ij}\}$ and $\{v_{ij}\}$.

Moreover the acceleration term $r_{ij}\cdot a_{ij} = r_{ij} \cdot(a_i - a_j)$ does involve the vector $a_i$ directly, but only the image of $a_i$ under the projection. Thus we could write the Newton Law as $$m_i a_i = F_{\text{effective}}+ \sum_j \frac{q_iq_j}{r_{ij}c^2} {\hat}}_{ij} (}}_{ij}}\cdot a_i), $$ and where $F_{\text{effective}}$ is the effective force. (We're being somewhat vague here).

These are all the terms that enter into the calculation of Weber's force law. It's more complicated than Newton's Gravitational equations, which depends only on the relative positions $\{r_{ij}\}$. Notice that $r_{ij}$, $r'_{ij}$, $r''_{ij}$ are scalar-valued, and the terms $r_i$, $v_i$, $v_{ij}$, and $a_{ij}$ are vector-valued.

Now the auxiliary function $f_i$ is ready-made for python's odeint. But the second right hand term is more interesting, as it involves the relative acceleration term directly, where $a_{ij}:=a_i - a_j$ by definition. We are looking to solve for $a_i$, for all $i$. The challenge is to:

  1. decouple the equations into independant ODEs
  2. try to rearrange the terms, see what symbolic manipulation can yield. This means writing it out in full (or in Wolfram Mathematica) and staring at it for a good while. Time for experimentation with the symbols.

Remark. An initial approximation might involve ignoring this second term. In fact we should rewrite our python code to reflect this direct factorization. Note that $r'_{ij}$ is still present in our formula, but we could directly use formula $r'_{ij}=\hat{\bf{r}}_{ij} \cdot v_{ij}$.