我正在模拟一个围绕原子核运行的电子。当然，带电粒子会向外辐射能量，因此它会撞击原子核。我的方法是评估库仑力并添加反作用力（由亚伯拉罕-洛伦兹公式给出）。运动方程为

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

在哪里$\beta = \frac{2}{3c^3}$在我的单位中，m = 1。重新排列这个def eq来解决是

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

我乘以$\frac{x}{r}$要么$\frac{y}{r}$分别获得 x 或 y 方向。实现如下

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()

我之前使用 odeint 并收到了错误消息Excess work done on this call，我知道可以使用可以处理事件的求解器来解决该错误消息，但是这些事件似乎只会在求解器中断之前终止它，但我希望它继续正常运行如果这有任何意义，就不要中断。使用 odeint 我得到以下情节

至少它在它应该做的原点完成，但它应该先绕过它几次。我想我的问题很明显，我怎样才能在不破坏的情况下解决颂歌？提前致谢。

编辑：我发现了问题（我认为）。当我们靠近中心时，1/r^3 变得太大并破坏了一切。但是我希望这一路走下去，我该如何解决这个问题？