我正在尝试解决一组耦合的非线性 ODE。唯一的因变量是一维空间坐标，我们称之为。现在，我已经设法近似消除了方程之间的一些耦合，这样第一个方程根本不依赖于第二个方程。我用 FEM 轻松解决了这个问题：满足边界条件，残差看起来不错，物理上有意义的结果。$x$

所以现在我只剩下一个非线性的二阶 ODE，其中非线性系数函数完全由 FEM 解决方案确定，例如

其中$A$、$B$、$C$和$D$是已知的，至少在 FEM 网格的分辨率下是已知的。

问题是，使用与以前相同的基本代码，解决方案似乎完全违反了边界条件，爆炸为$x \rightarrow 0$而不是接近零（甚至任何有限值，我似乎无关紧要选择）。认为这是我的代码的问题，我尝试了更精细的网格和有限差分而不是有限元，但没有成功。它每次都找到相同的荒谬结果。

最后，我将系数的数值数据放入 Mathematica 并尝试解决上述 BVP。在放弃之前，我让它运行了一段时间（可能是 15-30 分钟）。它从未找到解决方案。

所以最后，我的问题是：这个方程有可能根本无法用这些边界条件求解吗？例如，我尝试更改并且 Mathematica 几乎立即找到了解决方案，但它仍然违反了新的边界条件；不知何故，这并没有引发错误，但这让我怀疑这里可能会发生“更大”的事情。$f(1) = 0 \rightarrow f'(0) = 0$

我正在寻找为什么数值方法在这些情况下可能会失败以及如何缓解这些问题背后的一些数学直觉。

更新

这个方程是线性的（在这里更正以便评论仍然有意义）。不过，我认为这不是我问题的根源；我尝试使用具有适当刚度矩阵和载荷向量的线性求解器来求解，这应该很容易；这里，是第一个方程的 FEM 解。该解决方案满足边界条件，但“锯齿状”（想想覆盖在平滑函数上的锯齿图案）。现在我认为将第一个方程的 FEM 解代入第二个方程可能有问题，但我不知道为什么。$-f''(x) = g(x)$$g(x)$

还没有放弃。

更新 2

我想我可以安全地排除这段代码的问题；我正在使用 SymPy 生成源文件，因此不应该存在与更改我想要求解的方程的函数形式等相关的人为错误。作为我现在看到的示例，请参见下图。我对此的唯一直觉是代码“试图”满足边界条件（它们是齐次狄利克雷），但无论出于何种原因，解决方案基本上都想炸毁。所以我们最终得到了这种跨越两种可能性的“妥协”解决方案。$-\infty$

这甚至是远程正确的吗？在这一点上我改变什么似乎并不重要，细化网格只会增加锯齿形式的频率。