小线性系统的数值稳定显式解

计算科学 线性代数
2021-12-08 04:09:13

我有一个非齐次线性系统

Ax=b

在哪里A是一个真实的n×n矩阵与n4. 零空间A保证是零维的,所以方程有一个唯一的逆x=A1b. 由于结果进入 ODE 的右侧,我打算使用自适应方法求解,因此对于Ab. 由于这个要求和小维度,我想实现明确的公式A1b. 元素可以完全为零或采用非常不同的值。我的问题是这对您是否有意义,以及是否有已知的稳定表达式。我正在为 x86 系统编写 C 语言。

3个回答

在实现明确的公式之前,我会问自己一个问题:“值得吗?”:

  • 是否值得花时间编写、调试和验证这些显式公式,同时您可以轻松链接到使用经典高斯消元法的 BLAS+LAPACK?
  • 您期望获得稳定性吗?相反,我不认为编程显式公式(如克莱默规则)会给你更好的稳定性。
  • 您期望获得速度吗?您是否已经对整个程序进行了概要分析?解决这些 4x4 系统花费了多少时间?
  • 在一年的时间里,你改进你的模型并且你需要 5 个方程而不是 4 个的概率是多少?

我的建议:首先使用 BLAS/LAPACK 组合,看看它是否有效,分析整个程序,让学生实现明确的公式(抱歉,这里有点讽刺)并比较速度和鲁棒性。

我所知道的唯一明确的逆结果是Cramer's Rule,它最近被证明可以在O(n3)时间(如高斯消元法;但不确定主导因子前面的常数)。

的矩阵逆A是一个平滑函数A只要det(A)0,并且解肯定是的平滑函数,因此只要 ODE 的右侧是的平滑函数并且您避免秩不足的情况,我认为您的右手侧面会很光滑。(在这里,我用 smooth 表示“至少两次连续可微”。)xbxA

为了安全起见,最好确保在数值上也没有秩不足(即没有小的奇异值)。A

克莱默法则的问题在于它的稳定性属性是未知的,除了(它是前向稳定的,但不是后向稳定的)。(参见N. Higham的《数值算法的准确性和稳定性》,第 2 版。)它不被认为是可靠的算法;部分旋转的高斯消除 (GEPP) 受到青睐。n=2

我预计在 ODE 求解中使用 BLAS+LAPACK 执行 GEPP 的问题将是在隐式 ODE 方法中使用的任何有限差分。我知道人们已经解决了线性程序作为右手边评估的一部分,并且因为他们这样做很天真(只是将线性程序解决方案插入右手边,调用单纯形算法),他们大大降低了他们的准确性计算的解决方案,并大大增加了解决问题所需的时间。我的一个实验室伙伴想出了如何以更有效、更准确的方式解决这些问题。我得看看他的出版物是否已经发布。无论您选择使用 GEPP 还是 Cramer 规则,您都可能遇到类似的问题。

如果有任何方法可以为您的问题计算分析雅可比矩阵,您可能希望这样做以节省一些数字上的麻烦。与有限差分近似相比,它的评估成本更低,而且可能更准确。如果需要,可以在此处找到矩阵逆导数的表达式。评估矩阵逆的导数看起来至少需要两个或三个线性系统求解,但它们都将具有相同的矩阵和不同的右手边,因此它不会比单个线性系统贵很多解决。

如果有任何方法可以将计算出的解决方案与具有已知参数值的解决方案进行比较,我会这样做,以便您可以诊断是否遇到任何这些数值陷阱。

不确定这会有所帮助,但我只是认为当您谈论稳定解决方案时,您正在谈论近似方法。当您明确计算事物时,稳定性没有意义。也就是说,如果你想获得稳定性,你必须接受一个近似的解决方案。