solve_ivp - Overflow encountered in double_scalars
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 <PERSON>-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. 394
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.0c*3)/2.0)
def func(t, o): # [0x, 1y, 2vx, 3vy, 4ax, 5ay] r = np.hypot(o[0], o1) dodt = [o[2], o[3], o[4], o[5], Bo[0](-(1/r3.0) - o[4]/r), Bo1(-(1/r3.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], span1, 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 messageExcess 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? 509