Another post on Weber's force law, and not the last.
We return again to Weber force laws and two body problems, with special attention given to the possible attractive forces between identical negatively charged particles. This is one of the key predictions of Weber's force law, as addressed previously on this blog. The point is that a large amount of energy and work is required to push the particles through the critical radius r_c where suddenly the sign flips in the force term, or equivalently the object inherits a negative inertial mass. And this is the true meaning of E=mc^2 as we've explained.
But to progress this program forward we need to have ready integrated solutions to the equations of state. We also need to know that these integrated solutions are reasonably accurate. Practically we test their accuracy against the conservation of energy. Accurate solutions should maintain a constant total energy, namely $T+U$. Also one has conservation of angular momentum. That Weber's force law satisfies these conservation properties is proven by Weber himself repeatedly, and by Neumann, and again by AKT Assis in his recent publications.
We're returning to Weber two-body problem. Our earlier calculations were somewhat incorrect, since we did not account for the centre-of-mass frame properly. This has been corrected below.
It's curious somewhat that Weber's potential and force has caused much confusion and errors in physicists and mathematicians alike. It takes some care. Let's recall the basic derivation of the Weberian mass to borrow a phrase from Prof AKT Assis, especially in the following lecture slide.
Deriving the Critical Radius
We follow Assis' approach, as described in the above video. The idea is specific to the two-body system, where we have two particles $r_1$ and $r_2$ as an isolated system. The idea is to collect the $a_1$ terms, leaving the $a_2$ term within the force term. Thus Newton's 2nd Law applied to particle $r_1$ has the following form:
$$ \frac{q_1 q_2}} ~ {\hat{\bf{r}}}_{12} ~ \left(1-\frac{r'^2_{12}}{2c^2} - \frac{r_{12} a_2}{c^2} \right) = (m_1 - \frac{q_1 q_2}{r_{12} c^2}) ~a_1 .$$The left hand term is denoted $q_1 E$ in Assis' slide.
The critical radius $r_c$ is derived as the singular point, where we get a sign change, i.e. where $$r_c = \frac{q_1 q_2}{m_1 c^2}.$$
This formulation of the two-body problem is better suited to numerical integration.
Centre-Of-Mass Coordinates
We will be evaluating Weber's force law using Newton's 2nd Law in this centre-of-mass frame. This is linear algebra, but we made a mistake in our previous posts, so we here record the idea more formally.
We have $$R=\frac{1}{M} \sum m_j r_j$$ being the centre of mass in the initial Cartesian coordinate system, and also $$\frac{dR}{dt} = v_{CM} = \frac{1}{M} \sum m_j v_j.$$ The key observation is that the net force on $R$ is zero, i.e. $${\bf{F}} = M a_{cm} = 0.$$ [Verify this last statement]
Therefore $R(t) = v_{cm}.t$, and the idea of centre-of-mass coordinates is to "subtract" $v_cm$ from the reference frame. Thus we have a change of coordinates $$r'_i := r_i - R,$$ where $r_1', r_2', etc.,$ are going to be the position coordinates in the new centre-of-mass reference frame. The main observation is the trivial rearrangement $$r_i' = \frac{1}{M}\sum_{j} m_j r_{ij}.$$ This follows from the trivial $M-m_i = \sum_{j, ~j\neq i} m_j$.
Solving with GEKKO:
We would prefer to be integrating ODEs with odeint or more directly. However the form of Weber's force law implies that acceleration terms arise in the force, and we cannot use the usual tricks to isolate for each acceleration term.
Also since the computation of the Weber force requires relative accelerations between the particles, it has been necessary to record the accelerations in the state vectors. This is also somewhat inconvenient.
For example, say we being with initial state $s=[x,v,a]$. Newton's 2nd Law implies that the acceleration $a$ of the particle is equal to $1/m$ times the net force acting on the particle, i.e. total electric and non-electric forces. But we also find Weber's force $F$ depends on the relative accelerations (and distance and velocity) of the particles in the system. Therefore if only electrical forces are acting on the system, then what constraint does the acceleration term satisfy? The issue is that there's somewhat a circular definition at play which we need to carefully untie. Newton's law implies that the acceleration of the particle is equal to the net forces which act on the particle. However this net force itself depends on the relative acceleration of this particle with the system. So the question is: for a given time step $Δt$, how should the initial state $s$ propagate to $$s_{new} =s+Δt.(ds/dt)~~~?$$
By contrast, if we dealing with Newton's gravitational force with state $[x,v]$ (position and velocity) and law $F=ma$, then the propagation of any initial condition is more obviously $$[x_{new},v_{new}]=[x,v]+\Delta t.[v, \frac{1}{m}F].$$ This is the usual trick of setting $$dx/dt=v, ~~~dv/dt=a.$$ The above system can be directly implemented in odeint and numerically evaluated.
However in our case, with Weber's force law this is not immediately possible. Therefore we have resorted to GEKKO and its differential algebraic equation solver. But we do not like using other methods which offer limited visibility, because we cannot really be certain that the computation is correct!
!pip install gekko
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
c=1.0
def rho(x):
return np.dot(x,x)**(0.5)
# Component of Weber force which contains only position and velocities.
def f(state):
c=1.0
x, k = state
N=len(x)
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.
rhat=np.dot(1/rho(r), r) # rhat, unit vector.
Coulomb_factor = (+1)*k[i][-1]*k[j][-1]*(rho(r)**-2) # Coulomb factor, only term involving q values.
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.
s0 = (3/2)*((r_prime/c)**2)+(rho(dv)/c)**2
effective_Weber_factor = 1-s0
aux[i] = aux[i]+np.dot(Coulomb_factor*effective_Weber_factor, rhat)
return aux
def cm_frame(state):
x, k = state
N=len(x) # N-body system
N1=len(x[0]) # number of coordinates, expect N1 == 9.
output = np.zeros((N,N1))
output = np.array(output, dtype=object)
M=0
for i in range(N):
M = M + k[i][0]
for i in range(N):
for j in range(N):
if i==j:
pass
else:
dr = np.subtract(x[i], x[j])
output[i] = output[i] + np.dot((1/M)*k[j][0], dr)
return output, k
import numpy as np
from gekko import GEKKO
import matplotlib.pyplot as plt
c=1.0
# Weber two-body force.
def w2b(state):
c=1.0
x, k = state
N=len(x)
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 = (+1)*k[i][-1]*k[j][-1]*(rho**-2) # Coulomb factor, only term involving q values.
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
s0 = 1-(1/2)*(r_prime/c)**2
s1 = np.dot(r, x[j][6:9])*(c**-2)
s2 = k[i][0] - (k[i][1]*k[j][1])/(rho*c*c) ## this is the interesting term involving the critical radius.
factor_Weber = (s0 - s1)/s2 # scalar
aux[i] = aux[i]+np.dot(factor_Coulomb*factor_Weber, rhat)
return aux
def rho(x):
return np.dot(x,x)**(0.5)
# main function defining the effective Weber force.
def ef(state):
c=1.0
x, k = state
N=len(x)
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.
rhat=np.dot(1/rho(r), r) # rhat, unit vector.
Coulomb_factor = (+1)*k[i][-1]*k[j][-1]*(rho(r)**-2) # Coulomb factor, only term involving q values.
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.
s0 = (3/2)*((r_prime/c)**2)+(rho(dv)/c)**2
effective_Weber_factor = 1-s0
aux[i] = aux[i]+np.dot(Coulomb_factor*effective_Weber_factor, rhat)
return aux
# main function defining the Weber force.
def f(state):
c=1.0
x, k = state
N=len(x)
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 = (+1)*k[i][-1]*k[j][-1]*(rho**-2) # Coulomb factor, only term involving q values.
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
s0 = (1-(0.5)*(r_prime/c)**2)
s1 = (np.dot(dv, dv)- r_prime**2)*(c**-2) # scalar
s2 = np.dot(r, da)*(c**-2) # scalar
factor_Weber = s0 + s1 + s2 # scalar
aux[i] = aux[i]+np.dot(factor_Coulomb*factor_Weber, rhat)
return aux
# returns effective Weber mass. Scalar valued function.
def wm(state):
x,k=state
N=len(x)
for i in range(N):
k[i][0] - k[i]x[i]
def enw(state): ## here going to define the Effective Newton Weber Force
x,k=state
weber_force = f(state)
N=len(x)
return
def cm_frame(state):
x, k = state
N=len(x) # N-body system
N1=len(x[0]) # number of coordinates, expect N1 == 9.
output = np.zeros((N,N1))
output = np.array(output, dtype=object)
M=0
for i in range(N):
M = M + k[i][0]
for i in range(N):
for j in range(N):
if i==j:
pass
else:
dr = np.subtract(x[i], x[j])
output[i] = output[i] + np.dot((1/M)*k[j][0], dr)
return output, k
def gs(initial_state):
s = cm_frame(initial_state)
N=len(s[0])
m = GEKKO(remote=True)
m.time = np.linspace(0,5, 60) # time points
x = m.Array(m.Var,(N,9))
# initializing spatial values
for i in range(2):
for j in range(9):
x[i][j].value = s[0][i][j]
# setting constants.
c = np.zeros((2,2))
c = np.array(c, dtype=object)
for i in range(2):
for j in range(2):
c[i][0] = m.Const(s[1][i][0])
c[i][1] = m.Const(s[1][i][1])
# is eq0 correct? we see large accelerations, but the velocity is not increasing fast enough, so something seems wrong.
# eq0 = [ (c[i][0])*x[i][j+3].dt() == f([x, c])[i][j] for i in range(N) for j in range(3)]
eqA = [ x[i][j+3].dt() == w2b([x,c])[i][j] for i in range(N) for j in range(3)]
eqB = [ x[i][j].dt() == x[i][j+3] for i in range(N) for j in range(3) ]
#eqC = [ x[i][j].dt() == x[i][j+3] for i in range(N) for j in range(6) ]
m.Equation(eqA+eqB);
m.open_folder()
m.options.IMODE = 6
m.solve(disp=False)
return x
def IW(initial_state, repeat=1):
i = 0
input = initial_state
while i < repeat:
w=gs(input) # total solution integrated over time
terminal_state = [[w[i][j][-1] for j in range(9)] for i in range(2)] # this is terminal state of solution
# need to build the net_output array for the total integrated solution.
net_output = [] ## need initilize net_output to correct variable outside the scope of the if/else part.
net_output = np.asarray(net_output, dtype=object)
if i == 0:
net_output = w
else:
net_output = np.concatenate((net_output, w))
total_output = [terminal_state] + [input[1]] #
i = i+1
input = total_output
return net_output
def rprime(dx):
r = dx[0:3]
rho = np.dot(r,r)**(0.5) # rho = |r_ij|, scalar.
rhat = np.dot(1/rho, r)
dv = dx[3:6]
return np.dot(rhat, dv)
def U(dx):
r = dx[0:3]
rho = np.dot(r,r)**(0.5) # rho = |r_ij|, scalar.
return (rho**-1)*(1 - (c**-2)*(rprime(dx)**2))
def T(dx):
return rprime(dx)**2/2 ## missing a term?
#print(rprime([3.0,0.2,0, -1, 0, 0, 0, 0, 0]))
#print(U([3.0,0.2,0, -1, 0, 0, 0, 0, 0]))
#print(T([3.0,0.2,0, -1, 0, 0, 0, 0, 0]))
s0 = [[[0.2,0,0, -0.6,0.,0.1, 0,0,0], [0,0,0, 0,0,0.1, 0,0,0]], [[1.0, -1.0], [1.0, -1.0]] ];
a = IW(s0, 1);
#print(len(a[0][0]))
#print([a[0][0][j] for j in range(len(a[0][0]))])
#print(a[0][0][:])
#print(a[0])
#print(a[0][6][:])
s0 = [[[0.01,0,0, -0.6,0.,0.1, 0,0,0], [0,0,0, 0,0,0.1, 0,0,0]], [[1.0, -1.0], [1.0, -1.0]] ];
# w2b is renormalized vector representing the "effective force" term in Newton's 2nd Law.
#
print(
w2b(s0),"\n \n", f(cm_frame(s0))
)
# Okay, the sign of the force law _is_ correct.
# So we see a big difference between Weber's force, where the sign of the force is always repulsive for equal charge particles.
# Yet when we isolate the acceleration term and apply Newton's 2nd Law, then we find the "effective force" and the "effective mass"
# of the particle changes sign. To integrate the equations of motion, yes, we need to compute the "effective force". Need verify
# that our definition of w2b is correct.
fig = plt.figure(figsize=(20, 20))
ax = fig.gca(projection='3d')
plt.plot(a[0][0][:], a[0][1][:], a[0][2][:], "r", label='x0(t) with -1 charge')
plt.plot(a[1][0][:], a[1][1][:], a[1][2][:], "g", label='x1(t) with -1 charge')
#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) #? Does the rotation on perspective also rotate the axes!?
plt.show()
s1 = [[[1.2,0,0, 0,0,0, 0,0,0], [0.0,0,0, 0,0,0, 0,0,0]], [[1.0, -1.0], [1.0, -1.0]] ];
a = IW(s1, 1);
## oscillating position and velocity, but the velocity appears to be diverging.
## is this numerical error or correct integration of the equations of motion?
fig = plt.figure(figsize=(10, 10))
ax = fig.gca(projection='3d')
plt.quiver(a[0][0][:], a[0][1][:], a[0][2][:], a[0][3][:], a[0][4][:], a[0][5][:], color = 'r', label='x0(t) with -1 charge')
#plt.quiver(a[1][0][:], a[1][1][:], a[1][2][:], a[1][3][:], a[1][4][:], a[1][5][:], color = 'g', label='x1(t) with -1 charge')
#plt.quiver([a[0][1][0], a[0][1][1], a[0][1][2]], "g", label='x1(t) with -1 charge')
# 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) #? Does the rotation on perspective also rotate the axes!?
plt.show()
s0 = [[[1.01,0,0, 0,0,0, 0,0,0], [0.0,0,0, 0,0,0, 0,0,0]], [[1.0, -1.0], [1.0, -1.0]] ]
w2b(s0)
s0 = [[[0.9,0,0, 0,0.01,0, 0,0,0], [0.0,0,0, 0,-0.01,0, 0,0,0]], [[1.0, -1.0], [1.0, -1.0]] ];
a = IW(s0, 2);
# plottting
fig = plt.figure(figsize=(10, 10))
ax = fig.gca(projection='3d')
plt.plot(a[0][0][0], a[0][0][1], a[0][0][2], "r", label='x0(t) with -1 charge')
plt.plot(a[0][1][0], a[0][1][1], a[0][1][2], "g", label='x1(t) with -1 charge')
# 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) #? Does the rotation on perspective also rotate the axes!?
plt.show()
N=len(a[0][0][0])
print(N)
diff = [[a[0][0][j][k]-a[0][1][j][k] for j in range(9)] for k in range(N)]
print(
[dx[0] for dx in diff]
)
print(
[rprime(dx) for dx in diff]
)
# conservation of energy is not well managed. GEKKO is not doing a good job.
print(
[U(dx)+T(dx) for dx in diff]
)