首先得指出，你遇到的是刚性微分方程组的典型问题：当t增大时，a₀(t)和b₀(t)趋近于0，直接计算a(t)/a₀(t)这类项会触发浮点溢出，或者因为两个超大数值相减导致严重的精度丢失。不过好在你的系统有一个关键的守恒性质可以利用，再配合合适的求解策略，完全能解决稳定性问题。

一、先利用系统的守恒律简化方程 #

观察原方程组：

a'(t) = b(t)/b₀(t) - a(t)/a₀(t)
b'(t) = a(t)/a₀(t) - b(t)/b₀(t)

很容易发现a'(t) + b'(t) = 0，这意味着a(t) + b(t) = C，其中C是由初始条件确定的常数（比如初始时刻a(0)+b(0)的值）。

利用这个守恒量做变量替换，能直接把二维方程组降为一维，同时避免大项相减的精度问题：

• 令c(t) = a(t) - b(t)，那么a(t) = (C + c(t))/2，b(t) = (C - c(t))/2
• 代入原方程，推导得到c(t)的微分方程：

  c'(t) = (C - c(t))/b₀(t) - (C + c(t))/a₀(t)
  

这个一维方程的优势在于，它不再包含两个相近大值相减的项（原方程中a'(t)是两个超大数相减，极易丢失有效数字），计算精度和稳定性会大幅提升。

二、缩放变量避免除以极小值 #

如果不想降维，也可以通过变量缩放来规避直接除以极小的a₀(t)或b₀(t)：

• 定义新变量u(t) = a(t)/a₀(t)，v(t) = b(t)/b₀(t)，原方程的导数关系变为a'(t) = v(t) - u(t)，b'(t) = u(t) - v(t)
• 对a(t) = u(t)*a₀(t)求导，得到a'(t) = u'(t)*a₀(t) + u(t)*a₀'(t)，同理b'(t) = v'(t)*b₀(t) + v(t)*b₀'(t)
• 联立后得到关于u和v的方程：

  u'(t)*a₀(t) + u(t)*a₀'(t) = v(t) - u(t)
  v'(t)*b₀(t) + v(t)*b₀'(t) = u(t) - v(t)
  

这里的关键是，我们把极小的a₀(t)和b₀(t)移到了等式左边作为乘数，而不是除数，彻底避免了浮点溢出。另外，a₀'(t)和b₀'(t)有很简洁的解析式：

a₀'(t) = -20*a₀(t)*(1 - a₀(t))
b₀'(t) = -20*b₀(t)*(1 - b₀(t))

计算起来完全没有负担。

三、求解器选择：放弃显式RK，用刚性专用求解器 #

标准的显式Runge-Kutta方法（比如RK4）在处理刚性系统时效率极低，甚至会不稳定——因为你的系统在t较大时，特征值的绝对值会变得极大，显式方法需要极小的步长才能满足稳定性条件，计算成本会飙升。

更适合的是隐式求解器：

• BDF方法（向后差分公式）：这是工业界处理刚性ODE的首选，比如Python的SciPy中solve_ivp的method='BDF'选项，或者Matlab的ode15s，它能自适应调整步长，在保证稳定性的同时大幅减少计算量。
• 隐式Runge-Kutta方法：比如Radau IIA或Gauss方法，适合对精度要求较高的刚性问题，尤其是当系统非线性较强时，稳定性表现比BDF更好。

如果一定要用RK类方法，也可以选择半隐式RK，但整体来说，专用的刚性求解器是更优的选择。

四、额外的浮点精度优化技巧 #

• 确保使用双精度浮点数进行计算（大多数编程语言默认支持，但如果你的代码用了单精度，一定要切换），双精度的64位能更好地保留极小值附近的有效数字。
• 当t远大于0.6或0.3时，可以用近似公式计算a₀(t)和b₀(t)：比如当20*(t-0.6) > 10时，a₀(t) ≈ exp(-20*(t-0.6))，此时1/a₀(t) ≈ exp(20*(t-0.6))，这样能避免计算1/(1+exp(...))时的数值下溢，不过要注意只在t足够大时使用这个近似，避免过渡区域的误差。

内容的提问来源于stack exchange，提问作者JP-Ellis