I'm modeling an electron that orbits the nucleus. Of course, charged particles radiate away there energy so it'll crash into the nucleus. My approach has been to to evaluate the coulomb force and add the reaction force (given by the Abraham-Lorentz formula). The equation of motion would be

$F = ma = \frac{1}{r^2} - \beta\dot{a}$

where $\beta = \frac{2}{3c^3}$ in my units, and m =1. Rearranging this the def eq to solve is

$F = \dot{a} = (\frac{1}{r^2} - a) \frac{1}{\beta}$

and I multiply by $\frac{x}{r}$ or $\frac{y}{r}$ to get the x or y directions respectively. The implementation is below

import numpy as np
from scipy.integrate import odeint, solve_ivp
import matplotlib.pyplot as plt

c = 137.0
B = ((3.0*c**3)/2.0)

def func(t, o):  # [0x, 1y, 2vx, 3vy, 4ax, 5ay]
    r = np.hypot(o[0], o[1])
    dodt = [o[2], o[3], o[4], o[5], B*o[0]*(-(1/r**3.0) - o[4]/r), B*o[1]*(-(1/r**3.0) - o[5]/r)]
    return dodt

initial = [br, 0.0, 0.0, 1.0, -1.0, 0.0]  # [x, y, vx, vy, ax, ay]
span = [0.0, 10.0]
time = np.linspace(span[0], span[1], 10000)

# %% Solve differential equation
sol = solve_ivp(fun=func, t_span=span, t_eval=time, y0=initial, dense_output=True, max_step=0.05).sol
xe, ye, xve, yve, axe, aye = sol(time)

print(xe[-1], ye[-1])
# %% Plot states
plt.plot(xe, ye)
plt.plot(xe[-1], ye[-1], 'ro')
plt.plot(xe[0], ye[0], 'yo')
axx = plt.gca()

I was using odeint before and got the error message Excess work done on this call which from this post I know can be solved using a solver that can handle events, however it seems that these events only terminate the solver before it breaks, but I want it to continue as it shouldn't break to begin with if that makes any sense. Using odeint I get the following plot

enter image description here

At least it finishes at the origin which it's what it's supposed to do, but it should go around it a few times first. I guess my question is evident, how can I solve the ode without it breaking? Thanks in advanced.

Edit: I have found the problem (I think). When we get close to the center 1/r^3 becomes too large and breaks everything. However I would like this to go all the way, how do I get around this issue?


1 Answer 1


Are you sure about this expression:

B*o[0]*(-(1/r**3.0) - o[4]/r), B*o[1]*(-(1/r**3.0) - o[5]/r)

Doesn't the radiation force acting on the electron equal $$F_{rad} = \beta \, \frac{da}{dt}$$

The equations of motion, assuming they are in 2D which means we work with two component vectors $r, \, v, \, a \,\in \mathbb{R}^2$, are \begin{align} &\frac{dr}{dt} = v\\ &\frac{dv}{dt} = a\\ &\frac{da}{dt} = \left(\frac{B}{\,\,\|r\|^3}\right)\, r \, + \, B\,a \end{align} where $B = \frac{1}{\beta} = \frac{3c^3}{2}$.

Your interpretation is \begin{align} &\frac{dr}{dt} = v\\ &\frac{dv}{dt} = a\\ &\frac{da}{dt} = -\,B\left(\frac{r}{\|r\|^3} - \frac{a*r}{\|r\|}\right) \end{align}

  • $\begingroup$ Hi. Thanks for your reply. You are correct, I messed that up. But also, I've found that you can't really use the Abraham-Lorentz force because "it only tells you the time average of the parallel component—the average, moreover, over very special time intervals." That's from Griffiths. I at least it doesn't seem to me like you can. $\endgroup$
    – AJ09
    Commented May 5, 2020 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.