我一直在尝试解决与结构力学相关的非线性 PDE（准确地说是非线性 Timoshenko 梁）。我正在做h -refinement 和p -refinement 以达到解决方案。非线性方程组使用 Newton-Raphson 方法求解，我将我的解决方案（指定载荷下的横向位移）与基准解决方案（基准解决方案是使用高度p精化模型生成的，并且该基准解决方案被插值到解决方案使用样条插值的节点）。使用 -norm 和 -norm 计算相对误差。
。
代码是用 MATLAB 编写的。
说了这么多，我的问题是，解决方案的收敛（与基准）在一个限制后停止并开始围绕该限制振荡（参见附图以供参考-LFEM 是具有线性形状函数的 FEM，QFEM 是具有二次形状的 FEM函数，PFEM 是不同的p细化有限元）。我已经使用商业 FEM 包验证了该解决方案，并且它与它完美匹配。我已经梳理了整个代码的错误但无济于事。
错误的来源可能是什么？我观察到解决方案在某些有效数字之前都是准确的，即使在改进之后也不会改变。行为与预期一致，直到相对误差达到，然后问题开始。
有限元解决方案的总体误差通常包含许多分量。例如，这些可能包括 (i) 由有限网格和/或有限多项式次数引起的离散化误差，(ii) 求解定义有限元解的节点值的线性系统时产生的误差，(iii) 如果您的问题是非线性的，仅在有限多次迭代后终止非线性迭代（例如牛顿方案）所产生的误差，（iv）数值舍入所产生的误差；可能还有其他人。
为了使您的错误真正归零，您需要确保这些错误贡献中的每一个都归零。换句话说，当您细化网格时，您还需要更准确地求解线性系统，运行更多非线性迭代，并可能从单精度切换到双精度。如果您不控制这些错误源中的每一个，最终它们中的一个将支配所有其他错误源。例如，如果您只进行 3 次牛顿迭代，那么您求解每一次迭代的准确度并不重要：总体误差将取决于您在 3 次牛顿迭代后停止的事实。
最后，您的“参考”解决方案也有可能不完全准确。例如，如果您从一个文件中读取它，该文件中存储了 6 位精度的参考解的节点值，您不能期望您的数值解和这个“参考”解之间的误差小于仅仅是因为参考解决方案的最后一位不再准确。
对于你的问题，我真的没有明确的答案。但是，经过一些实验，我相信您关于刚度矩阵条件数的结论是正确的。
我也相信您可以使用迭代求解器来提高准确性。一个非常好的预处理器是使用直接解中的因子。在 Octave/MATLAB 中，您可以使用以下代码：
lt = chol(K);
u = pcg(K, rhs, 1e-16, 10, lt', lt);
这基本上是一种迭代细化的形式。
这不是这个问题的最终结论，但它为更多实验指明了方向。