计算科学 - 应用有限差分方案时矩阵向量乘积的精度损失 - 吾爱随笔录

应用有限差分方案时矩阵向量乘积的精度损失

计算科学线性代数有限差分数值分析精确准确性

2021-12-23 03:46:51

我正在应用 6 阶有限差分微分方案，如http://www.scholarpedia.org/article/Method_of_lines/example_implementation/dss006中所示

该方案似乎存在严重的数字/浮点问题。结果对与 FD 矩阵预乘的向量中的值非常敏感。

例如，具有相同值的 10*1 向量的梯度全为零。例如，应用下面的 10*10 矩阵

   -1764        4320       -5400        4800       -2700         864        -120           0           0           0
    -120        -924        1800       -1200         600        -180          24           0           0           0
      24        -288        -420         960        -360          96         -12           0           0           0
     -12         108        -540           0         540        -108          12           0           0           0
       0         -12         108        -540           0         540        -108          12           0           0
       0           0         -12         108        -540           0         540        -108          12           0
       0           0           0         -12         108        -540           0         540        -108          12
       0           0           0          12         -96         360        -960         420         288         -24
       0           0           0         -24         180        -600        1200       -1800         924         120
       0           0           0         120        -864        2700       -4800        5400       -4320        1764

到 10x1 向量25545.007*ones(10,1)

结果 [6.1118e-10, -2.3574e-09, -1.266e-09, -3.638e-10, -3.638e-10, -3.638e-10, -3.638e-10, 2.0082e-09, - 6.1118e-10，-9.8225e-09]

而正确答案是 10x1 零向量。同一组向量上的二阶导数矩阵产生的结果更糟，其数量为 1e3 而不是零！

当我在一个稍微不同的向量上尝试这个相同的一阶导数 FD 矩阵时25545.007*ones(10,1)，它会起作用并产生正确的答案（即零向量）。这些结果已在 MATLAB、Python(Numpy) 和 Mathematica 上使用标准双精度算术进行了双重检查。每个软件的确切数值不同，但数量级与上述结果相似。

为什么该方案对仅相差 1e-3 的向量如此敏感？我真的需要高精度，这样的数值问题意味着我在下游计算中的解决方案值变得不可靠，并且对整体解决方案没有信心。

有没有更好的方法来实现这一点？对于这种情况，也许存在替代的推荐方案？由于项目的复杂性及其固有的缓慢性，我无法引入可变精度算术库。

2个回答

如果我理解正确，您使用的是 LLNL SUNDIALS 套件中的IDA求解器。

如果您可以访问计算 DAE 系统残差的源代码，那么您可以考虑自动微分和IDAS，这基本上是具有前向和伴随灵敏度分析功能的 IDA。

有一个用户手册和一个完整的示例列表。

如果你想坚持数值微分，我认为这个 na.stackexchange 帖子是一个很好的起点。

这可能只是一个舍入错误问题。微分矩阵的每一行总和必须为零。减去向量的平均值可以得到更好的结果

x = 25545.007*numpy.ones((10,1))
A.dot(x)
array([[  0.00000000e+00],
       [ -7.45058060e-09],
       [  0.00000000e+00],
       [ -1.01863407e-10],
       [ -1.01863407e-10],
       [ -8.14907253e-10],
       [ -9.31322575e-10],
       [ -3.72529030e-09],
       [  3.72529030e-09],
       [  1.49011612e-08]])
xa = x.mean()
y = x - xa
A.dot(y)
array([[ 0.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 0.]])

另一种选择是这个。如果导数由下式给出

v = D u, v_{i} = \sum_{j = 1}^{N} D_{i j} u_{j}

$v = Du, \qquad v_i = \sum_{j=1}^N D_{ij} u_j$ 然后计算为

v_{i} = \sum_{j = 1}^{N} D_{i j} (u_{j} - u_{i}), 1 \leq i \leq N

$v_i = \sum_{j=1}^N D_{ij} (u_j - u_i), \qquad 1 \le i \le N$ 我希望这可以减少舍入误差。

其它你可能感兴趣的问题

上一篇“具有对角线优势”的稠密矩阵的预处理器下一篇数值查找格林函数