Newton-Krylov 何时不是合适的求解器?

计算科学 Python 非线性方程 泊松 scipy
2021-12-09 22:39:52

最近,我一直在比较scipy 中的不同非线性求解器,并且对 Scipy Cookbook 中的 Newton-Krylov 示例印象特别深刻,他们在大约 20 行代码中求解了具有非线性反应项的二阶微分方程。

我修改了示例代码来求解半导体异质结构的非线性泊松方程(也称为泊松-玻尔兹曼方程,参见这些注释中的第 17 页),其形式为,

d2ϕdx2k(x)(p(x,ϕ)n(x,ϕ)+N+(x))=0

(这是传递给求解器的残差函数。)

这是一个静电问题,其中n(x,ϕ)p(x,ϕ)是形式的非线性函数ni(x)e(Ei(x,ϕ)Ef). 这里的细节并不重要,但重点是非线性函数随ϕ所以残差函数可以在很大的范围内变化(1061016)略有变化ϕ.

我用 scipy 的 Newton-Krylov 数值求解这个方程,但它永远不会收敛(实际上它总是会在计算雅可比行列时报告错误)。我从Newton-Krylov求解器切换到fsolve(基于 MINPACK hybrd),它第一次工作!

Newton-Krylov 不能很好地解决某些问题是否有一般原因?输入方程是否需要以某种方式进行调节?

也许需要更多信息才能发表评论,但您为什么认为 fsolve 在这种情况下有效?

1个回答

您可能会遇到两个问题。

病态

首先,这个问题是病态的,但如果你只提供一个残差,Newton-Krylov 会通过对残差进行有限差分来丢弃一半的有效数字,以获得雅可比行列式的作用:

J[x]yF(x+ϵy)F(x)ϵ

如果您提供解析雅可比行列式,则可以保留所有数字(例如,双精度 16)。Krylov 方法也依赖于内积,所以如果你的雅可比行列式不适用于1016,它实际上是奇异的,并且 Krylov 可以停滞或返回错误的解决方案。这也可以防止直接求解器收敛。有时您可以使用多重网格方法将其均质化为具有易处理条件的粗网格。当一个问题不能用更好的条件来表述时,它可能值得在四精度下工作。(例如,这由 PETSc 支持。)

请注意,同样的问题也适用于准牛顿方法,尽管没有有限差分。对于非紧算子问题(例如微分方程)的所有可扩展方法都必须使用雅可比信息进行预处理。

很可能fsolve要么不使用雅可比信息,要么使用狗腿法或转移以通过“梯度下降”方法取得进展,尽管本质上是奇异的雅可比(即,有限差分会从有限精度算术)。这是不可扩展的,并且fsolve随着问题大小的增加可能会变慢。

全球化

如果线性问题得到准确解决,我们可以排除与线性问题(Krylov)相关的问题,而专注于非线性问题。局部最小值和非平滑特征收敛缓慢或导致停滞。Poisson-Boltzmann 是一个平滑模型,所以如果你从一个足够好的初始猜测开始,牛顿将二次收敛。大多数全球化策略都涉及某种延续,以便为最终迭代产生高质量的初始猜测。示例包括网格延续(例如,完全多重网格)、参数延续和伪瞬态延续。后者通常适用于稳态问题,并提供了一些全局收敛理论,参见Coffey、Kelley 和 Keyes (2003)搜索出现了这篇论文,这可能对你有用:Shestakov、Milovich 和 Noy (2002)使用伪瞬态延拓和有限元法求解非线性 Poisson-Boltzmann 方程伪瞬态延续与 Levenberg-Marquardt 算法密切相关。

进一步阅读