How to obtain and form a 1st order differential equation for leapfrog integration from second order one in this example of coulomb drag

Question

I am currently doing a computational physics homework which asked us to use leapfrog to give the relations between timevelocities and time-distance of these two objects.

The full question is as follows:

Consider two parallel wire with distance d. A negative charged object(m1) is moving along the wire with initial velocity v0. One positive charged object(m2) is located on another wire with zero initial velocity and located at the original point(x = 0). The initial distance between two objects(x), velocity(v0), charges and(e1,e2) masses(m1,m2) are arbitrary. Relativity effect can be ignored.

1. Write down the equations of motion for these two objects for this problem and convert it from a second-order equation to two-rst order ones. Write down a program using the time-reversal symmetric methods such as leapfrog or Verlet method to give the relations between timevelocities and time-distance of these two objects.(please use initial conditions: m1 = 1; m2 = 1; e1 = −1; e2 = 1; d = 1; x = 100; v0 = −5)
2. Keep the original distance, masses and charges unchanged(m1 = 1; m2 = 1; e1 = −1; e2 = 1; d = 1; x = 100), please plot the nal velocity of m2 with dierent v0 (v0 ∈ (−10,0)). (here, you may need to enlarge the total time to nd the stable nal velocity of m2, roughly should be δv < 10−4). In the distance gure, we can always nd that m2 rst move along positive direction and then move along negative direction. Please nd the maxinum distance(D2max) of m2 in the positive direction at dierent v0 and plot this relation between v0 and D2max

I have already found the two 2nd order differential equation:

And by substituting the initial value as required, we obtained the following result:

And I have already convert the two equations into first order:

First of all, I would like to ask whether or not my steps for now are correct or not, since I am really skeptical on whether I can just directly convert it to first order one since the upsidedown L in the denominator still contains the variable x_1 and x_2

Secondly, since there are variable x_1 and x_2 in the upsidedown L, how can I change the equation to the form that I can use in the modified function of rk2a() in http://www.math-cs.gordon.edu/courses/mat342/python/diffeq.py

I have modified the code in rk2a() function to the leapfrog form:

def leapfrog(f, x0, t):
n = len( t )
x = np.array( [ x0 ] * n )
for i in range( n - 1 ):
    h = t[i+1] - t[i]
    k1 = h * f( x[i], t[i] ) / 2.0
    if (i%2 == 0):
        x[i+1] = x[i] + h * f( x[i] + k1, t[i] + h / 2.0 )
    else:
        k2 = h * f( x[i], t[i] )
        x[i+1] = k1 + h * f( x[i] + k2, t[i] + h / 2.0 )

return x

Thirdly, I would like to ask whether my modification is correct or not(Hopefully it should be correct i guess, but I don't have any equation to test it with since above, I have trouble changing it to the form required by the function)

And the modification is based on the equatuion by my lecture notes:

Answer 1

You have already found the system of two differential equations (I am assuming that by $\Gamma_{12}$ you mean $\Gamma_{12} = \sqrt{(x_1-x_2)^2 + d^2\,}$ ) which is \begin{align} &m_1 \frac{d^2 x_1}{dt^2} = \frac{e_1e_2}{4\pi \varepsilon_0} \, \frac{x_1-x_2}{\big(\,(x_2-x_1)^2 + d^2\,\big)^{\frac{3}{2}}}\\ &m_2 \frac{d^2 x_2}{dt^2} = \frac{e_1e_2}{4\pi \varepsilon_0} \, \frac{x_2-x_1}{\big(\,(x_1-x_2)^2 + d^2\,\big)^{\frac{3}{2}}} \end{align} According to the assumption of the problem $m_1=m_2=1, \,\, e_1=-1,\,\,e_2=1$ and $d=1$. Consequently \begin{align} &\frac{d^2 x_1}{dt^2} = k_0 \, \frac{x_2-x_1}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}}\\ &\frac{d^2 x_2}{dt^2} = k_0 \, \frac{x_1-x_2}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}} \end{align} where $k_0 = \frac{1}{4\pi \varepsilon_0}$ (seems like you re missing this factor). Furthermore, this system reduces from two equations to one, due to the translation invariance (i.e. there is a one parameter group of translation symmetries along the $x$ axis) which provides a conserved quantity, by Noether's theorem. But we can carry out the reduction by hand, without too much fancy noise and big words, as follows: subtract the second equation from the first one and obtain $$\frac{d^2 x_1}{dt^2} - \frac{d^2 x_2}{dt^2} = - k_0 \, \frac{x_1-x_2}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}} - k_0 \, \frac{x_1-x_2}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}}$$ which turns into $$\frac{d^2}{dt^2}(x_1 - x_2) = \frac{- 2k_0 \,(x_1-x_2)}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}}$$ Next, set $x = x_1-x_2$ $$\frac{d^2 x}{dt^2} = - \, \frac{2k_0 \, x}{\big(\,1 + x^2\,\big)^{\frac{3}{2}}}$$ The latter equation is the differential equation that describes the time-evolution of the relative position of the two particles. I believe this is the equation that is mentioned in part 1 of your questions (notice is says "and convert it from a second-order equation to two-rst order ones).

Furthermore, if you sum the two equations, you obtain the relation $$\frac{d^2 x_1}{dt^2} + \frac{d^2 x_2}{dt^2} = - k_0 \, \frac{x_1-x_2}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}} + k_0 \, \frac{x_1-x_2}{\big(\,(x_1-x_2)^2 + 1\,\big)^{\frac{3}{2}}} = 0$$ i.e. $$\frac{d^2}{dt^2}(x_1 + x_2) = 0$$ and so by integrating two times consequtively $$\frac{d x_1}{dt}(t) + \frac{dx_2}{dt}(t) = v_0 + 0 = -5$$ $$x_1(t) + x_2(t) = x_1(0) + x_2(0) - 5\,t = 100 + 0 - 5\,t = 100 - 5\,t$$ Combining the latter equation with the equation $x_1(t) - x_2(t) = x(t)$ allows you to express everything in terms of $x = x(t)$ \begin{align} &x_1(t) = \frac{1}{2}\left(100 - 5\, t + x(t)\right)\\ &x_2(t) = \frac{1}{2}\left(100 - 5\, t - x(t)\right)\\ &\\ &v_1(t) = \frac{dx_1}{dt}(t) = \frac{1}{2}\left(- 5 + v(t)\right)\\ &v_2(t) = \frac{dx_2}{dt}(t) = \frac{1}{2}\left(- 5 - v(t)\right)\\ \end{align} where $v(t) = \frac{dx}{dt}(t)$. On a side note, the differential equation for $x$ has a conserved energy integral, so it is in fact completely integrable (solvable in possibly implicit functions and integrals). Anyway, the second order differential equation for $x$ turns into a system of two first order equations \begin{align} &\frac{dx}{dt} = v\\ &\frac{dv}{dt} = - \,\frac{ 2k_0 \, x\, }{\big(\,1 + x^2\,\big)^{\frac{3}{2}}} \end{align} The way I know leapfrog, one looks at the system as a sum of two exactly solvable summands, i.e. \begin{align} \frac{d}{dt}\begin{bmatrix}x\\v\end{bmatrix} = \begin{bmatrix}v\\0\end{bmatrix} + \begin{bmatrix} 0\\ - \, \frac{ 2k_0 \, x\, }{\big(\,1 + x^2\,\big)^{\frac{3}{2}}}\end{bmatrix} \end{align} and based on this decomposition, split it into two exactly solvable systems, i.e. \begin{align} &\frac{dx}{dt} = v\\ &\frac{dv}{dt} = 0 \end{align} and \begin{align} &\frac{dx}{dt} = 0\\ &\frac{dv}{dt} = - \,\frac{ 2k_0 \, x\, }{\big(\,1 + x^2\,\big)^{\frac{3}{2}}} \end{align} The solutions to the first system look like \begin{align} &x(t) = x_0 + v_0\,t\\ &v(t) = v_0 \end{align} and the solutions of the second system are \begin{align} &x(t) = x_0\\ &v(t) = v_0 - t\,\frac{ 2k_0 \, x_0\, }{\big(\,1 + x_0^2\,\big)^{\frac{3}{2}}} \end{align} Leapfrog is basically the carefully chosen composition of these two sets of solutions: within one time step $h$, starting from the current $x_i$ and $v_i$, follow the first solution for time $h/2$, then the second solution for time $h$ and finally again the first solution for time $h/2$, finally arriving at $x_{i+1}$ and $v_{i+1}$: \begin{align} &x_{i+1/3} = x_i + \frac{h}{2}\,v_i \\ &v_{i+1/3} = v_{i} \\ &\\ &\\ &x_{i+2/3} = x_{i+1/3}\\ &v_{i+2/3} = v_{i+1/3} - h\,\frac{ 2k_0 \, x_{i+1/3}\, }{\big(\,1 + x_{i+1/3}^2\,\big)^{\frac{3}{2}}} \\ &\\ &\\ &x_{i+1} = x_{i+2/3} + \frac{h}{2}\,v_{i+2/3}\\ &v_{i+1} = v_{i+2/3} \end{align}

So in Matlab code (I do not master python that much) you can write

X1 = zeros(1,n); %initialization of a one row array of length n
V1 = zeros(1,n); %initialization of a one row array of length n
X2 = zeros(1,n); %initialization of a one row array of length n
V2 = zeros(1,n); %initialization of a one row array of length n
x = 100; 
v=-5;
t=0;
X1(i) = 100;
X2(i) = 0 ;
V1(i) = -5;
V2(i) = 0;
for i = 2:n
   t = t + h; 
   x = x + v*h/2;
   v = v - 2*k0*h*x/(1+x^2)^(3/2);
   x = x + v*h/2;
   X1(i) = (100 - 5*t + x)/2;
   X2(i) = (100 - 5*t - x)/2;
   V1(i) = (- 5 - v)/2;
   V2(i) = (- 5 + v)/2;
end