如何使用Boost odeint数值求解含奇点的微分方程并避免除零错误？

针对柯西-欧拉方程、勒让德方程这类带正则奇点的微分方程，要避免除零错误并使用Boost odeint数值求解，可采用以下几种方法：

1. 变量替换消去奇点 #

通过变量替换将原方程转换为无奇点的常系数ODE，再直接用Boost odeint求解：

• 柯西-欧拉方程：对于方程 x²y'' + axy' + by = 0，令 t = ln(x)（仅适用于x>0），代入后可将原方程转换为常系数线性ODE：

  y''_t + (a-1)y'_t + by = 0
  

  设y₁=y，y₂=y'_t，即可得到标准一阶方程组，直接用Boost odeint的任意积分器求解，最后将t转换回x即可。
• 勒让德方程：对于方程 (1-x²)y'' - 2xy' + n(n+1)y = 0，令 x = cosθ，代入后可转换为关于θ的ODE：

  d²y/dθ² + cotθ dy/dθ + n(n+1)y = 0
  

  该方程的奇点性质被弱化，可通过Boost odeint的自适应步长积分器稳定求解。

2. 保留原方程形式，用自适应步长积分器处理 #

转换一阶方程组时不直接除以含奇点的系数，而是保留原式结构，借助自适应步长积分器自动规避数值问题：
以勒让德方程为例，设y₁=y，y₂=y'，原方程可整理为：

(1-x²)y₂' = 2x y₂ - n(n+1)y₁

在x≠±1时，y₂' = (2x y₂ - n(n+1)y₁)/(1-x²)。由于解析解在x=±1处有界，分子分母会同时趋近于0，此时使用Boost odeint的自适应步长积分器（如runge_kutta_cash_karp54），并设置合适的绝对误差abs_error和相对误差rel_error阈值，积分器会自动缩小奇点附近的步长，数值上可稳定计算0/0的极限值，避免除零错误。

3. 结合解析近似处理奇点区域 #

利用正则奇点处的解析解性质，在奇点附近切换为解析近似计算，远离奇点后再用数值积分：
比如柯西-欧拉方程在x=0附近，解析解为y = C₁x^r₁ + C₂x^r₂（r₁、r₂为特征方程根），当|x| < ε（ε为自定义小阈值）时，直接用该解析表达式计算y和y'；当|x| ≥ ε时，再切换回Boost odeint的数值积分流程，彻底避免奇点处的除法操作。

内容的提问来源于stack exchange，提问作者GregReese