隐式Runge-Kutta（Radau IIA）求解刚性方程不收敛问题求助

针对你遇到的Newton迭代发散/停滞、收敛效果不佳的问题，从刚性方程特性与隐式RK实现细节出发，给出以下排查方向：

1. Newton迭代初始猜测优化 #

当$\delta$极小时（如$1e-6$），区间$[0,2/\delta]$长达$2e6$，解会快速趋近于稳态$y=1$，后续几乎无变化。若沿用前一步解作为Newton迭代初始猜测，接近稳态时雅可比矩阵易接近奇异，直接导致迭代发散。

• 优化方案：解接近稳态时，改用当前步的$y$值作为常数初始猜测，或用前3步解做线性外推生成初始值，避免猜测值偏离真实解过多。

2. 雅可比矩阵的计算与数值稳定性 #

隐式RK的Newton迭代依赖雅可比矩阵$J=I - hAf_y$，其中$f_y$是$f(y)=y^2–y3$的导数（即$2y-3y^2$）。当$y\to1$时，$f_y=-1$，雅可比矩阵条件数会显著增大，步长$h$较大时数值稳定性急剧下降。

• 检查要点：
  • 是否在每个Newton迭代步都重新计算$f_y$？
  • 求解线性方程组时是否使用带选主元的LU分解？避免直接求逆导致的奇异性问题。
• 优化方案：接近稳态时，显式计算雅可比矩阵并启用选主元LU分解，必要时可给对角元素添加极小扰动（如$1e-12$），避免矩阵奇异。

3. 自适应步长控制逻辑调整 #

$\delta=1e-2$时区间长度为200，解从初始值（通常为$y(0)=\delta$）快速上升到1再进入稳态，若步长控制仅依赖局部误差估计，易出现两个极端：

• 快速变化阶段步长过大，Newton迭代无法收敛；
• 稳态阶段步长过小，浪费计算资源且易因机器精度限制导致迭代停滞。
• 检查要点：
  • Newton迭代发散时是否及时减半步长并重试？
  • 是否针对刚性问题设置了最大步长限制（如不超过$1/\delta$）？
• 优化方案：
  • 新增迭代失败检测逻辑，一旦Newton迭代停滞/发散，立即将步长减半后重新计算；
  • 解进入稳态后，固定步长或放宽步长上限，同时调整局部误差估计的阈值。

4. Newton迭代收敛阈值的自适应设置 #

固定的收敛阈值在刚性问题中易失效：阈值过严，稳态阶段因残差已接近机器精度，迭代会陷入停滞；阈值过松，快速变化阶段会引入过大局部误差。

• 检查要点：是否使用相对残差作为停止条件？比如判断norm(residual) < tol*norm(y_guess)，而非固定绝对阈值。
• 优化方案：采用自适应阈值，快速变化阶段用较小相对阈值（如$1e-8$），稳态阶段切换为较大相对阈值（如$1e-6$）。

5. 基础实现细节核查 #

最后确认核心代码的正确性：

• Radau IIA的系数矩阵$A$是否准确（以3阶为例，$A=\begin{bmatrix}5/12 & -1/12 \ 3/4 & 1/4\end{bmatrix}$，节点为$[1/3,1]$）；
• 隐式方程残差计算是否正确：残差应为$Y - y_{\text{prev}} - hAf(Y)$，其中$Y$是阶段值向量；
• Newton迭代更新公式是否正确：$Y = Y - J\backslash\text{residual}$，确保线性方程组求解逻辑无误。

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