In previous posts we have been focused on modelling N-body systems under Weber's force law. Concretely we have taken $N=2,3$. This post focusses on $N=2$. The initial challenge of these simulations was finding a suitable DAE (Differential Algebraic Equation) solver. This is because Newton's 2nd Law and Weber's force lead to implicit ordinary differential equations (nonlinear, second order). We have found GEKKO to be suitable for our purposes, being easy to install, well maintained by Prof. Dr. Scott Hedengren, and easy to use. Of course GEKKO has limitations based on the input size. In our case even two particles involves 18 variables and the ODE system consists of 12 equations. This is somewhat demanding for usual computers to handle.

!pip install gekko

Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/
Requirement already satisfied: gekko in /usr/local/lib/python3.7/dist-packages (1.0.5)
Requirement already satisfied: numpy>=1.8 in /usr/local/lib/python3.7/dist-packages (from gekko) (1.21.6)

import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt

# main function defining the Weber force.
def f(state):
  c=1.0
  N=len(state)
  x = state
  aux = np.zeros((N,3))
  aux = np.array(aux, dtype=object)
  for i in range(N):
    for j in range(N):
      if i==j:
        pass
      else:
        r=np.subtract(x[i][0:3], x[j][0:3])  # r_ij, vector.
        rho=np.dot(r,r)**(0.5) # rho = |r_ij|, scalar.
        rhat=np.dot(1/rho, r) # rhat, unit vector.
        factor_Coulomb = x[i][-1]*x[j][-1]*(rho**-2) # Coulomb factor
        dv = np.subtract(x[i][3:6], x[j][3:6]) # dv= v_{ij}, vector.
        r_prime=np.dot(rhat, dv) # r' = rhat \cdot dv, scalar.
        da = np.subtract(x[i][6:9], x[j][6:9]) #da = a_{ij}, vector
        s1 = (np.dot(dv, dv)- r_prime**2)*(c**-2) # scalar
        s2 = np.dot(r, da)*(c**-2) # scalar
        factor_Weber_1=(1-(0.5)*(r_prime/c)**2 + s1 + s2 ) # scalar

        aux[i]=aux[i]+np.dot(factor_Coulomb*factor_Weber_1, rhat)
       
  return aux #returns a force vector, i.e. three-vector.


# returns the centre of mass of the initial state. 
def cm(state):
  xx = state
  N=len(state)
  M=np.zeros(N)
  state_cm=[0,0,0,0,0,0,0,0,0]
  
  for i in range(N):
    M[i] = np.abs((xx[i][-1])**-1)
    state_cm= np.add(state_cm, np.dot(M[i], xx[i][0:9]))

  return np.dot(state_cm, (sum(M))**-1)

# better results obtained when remote=True
m = GEKKO(remote=True)  

# the following line needs be adjusted frequently 
# if m.solve() doesn't successfully terminate.
m.time = np.linspace(0,15, 250) # time points

#print("m.time returns: \n", m.time)

# create GEKKO variables
x = m.Array(m.Var,(2,10)) 

#

# initial positions
x[0][0].value = +0.2
x[1][0].value = -0.2

#x.value[:][0:3]=terminal_state



#x[2][2].value = +1.1

# initial velocities
x[0][3].value=-0.6
x[0][5].value=0.1
x[1][3].value=+0.7
x[1][5].value=-0.1


# defining the charges and inertial mass
x[0][-1].value=-1.0  # q0/m0
x[1][-1].value=-1.0  # q1/m1
#x[2][-1].value=+0.01  # q2/m2

print("The initial states are: \n", x, "\n")

# equations defining the ODE system
eq0 = [ x[i][j+3].dt() == f(x)[i][j] for i in range(2) for j in range(3)] 
eq1 = [ x[i][j].dt() == x[i][j+3] for i in range(2) for j in range(3)]

print(len(eq0+eq1))

m.Equation(eq0+eq1);

The initial states are: 
 [[0.2 0 0 -0.6 0 0.1 0 0 0 -1.0]
 [-0.2 0 0 0.7 0 -0.1 0 0 0 -1.0]] 

12

Run gekko solver:

m.open_folder()
m.options.IMODE = 6
m.solve(disp=False)

# the centre of mass position and velocity.
s_init=[]
N=2 # two body system.
for i in range(N):
  s_init=s_init + [[x[i][j][0] for j in range(10)]]

cm_s = cm(s_init)
v_cm = cm_s[3:6]
RR_cm = [np.add(np.dot(v_cm, time), cm_s[:3]   ) for time in  m.time ]
Rt_cm=np.transpose(RR_cm)

print("the initial state is:", s_init)  ## returns initial state from solution.
print("cm of initial state is:", cm_s) ## returns cm from initial state. This function is defined earlier.
print("velocity of cm is:", v_cm)  ## v_cm is constant with respect to time (verify with computation)

the initial state is: [[0.2, 0.0, 0.0, -0.6, 0.0, 0.1, 0.0, 0.0, 0.0, -1.0], [-0.2, 0.0, 0.0, 0.7, 0.0, -0.1, 0.0, 0.0, 0.0, -1.0]]
cm of initial state is: [0.   0.   0.   0.05 0.   0.   0.   0.   0.  ]
velocity of cm is: [0.05 0.   0.  ]

# and subtracting the centre of mass to prepare for the plotting.

x0_cm = np.subtract(x[0][0], Rt_cm[0])
y0_cm = np.subtract(x[0][1], Rt_cm[1])
z0_cm = np.subtract(x[0][2], Rt_cm[2])

x1_cm = np.subtract(x[1][0], Rt_cm[0])
y1_cm = np.subtract(x[1][1], Rt_cm[1])
z1_cm = np.subtract(x[1][2], Rt_cm[2])

l = len(x[0][0]);
vx0 = [x[0][3][j]-v_cm[0] for j in range(l) ]
vy0 = [x[0][4][j]-v_cm[1] for j in range(l) ]
vz0 = [x[0][5][j]-v_cm[2] for j in range(l) ]

vx1 = [x[1][3][j]-v_cm[0] for j in range(l) ]
vy1 = [x[1][4][j]-v_cm[1] for j in range(l) ]
vz1 = [x[1][5][j]-v_cm[2] for j in range(l) ]


# Creating new plot.
fig = plt.figure(figsize=(30, 30))
ax = fig.gca(projection='3d')

#ax0 = fig.add_subplot(111, projection='3d')
#plt.gca().patch.set_facecolor('black')

plt.plot(x0_cm, y0_cm, z0_cm, "r:", label='x0(t) with -1 charge')
plt.plot(x1_cm, y1_cm, z1_cm, "g:", label='x1(t) with -1 charge')
plt.plot(x[0][0], x[0][1], x[0][2], "r")
plt.plot(x[1][0], x[1][1], x[1][2], "g")

# Add Title
plt.title("Two body system of identical negatively charged particles interacting with respect to Weber's electrodynamical force law.") 
plt.legend(loc='best')
# Display
plt.axis('on')
ax.view_init(30, 20)
plt.show()