solve_ivp - Overflow encountered in double_scalars

Question

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()
axx.set_aspect('equal')
plt.show()

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

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?

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}