C和MATLAB中的矩阵乘法，不同的结果

Question

我正在使用 4Rungekutta 来解决 MATLAB 和 C 中的 DGL(Ax + Bu = x_dot)， A 是 5x5，x 是 5x1，B 5x1，u 1x1，u 是正弦函数的输出（2500 点）， MATLAB 和 C 中 4Rungekutta 的输出在第 45 次迭代之前都是相同的，但是在 4Rungekutta 的第 45 次（2500 次迭代中）迭代时，4Rungekutta 的第 2 步的 A*x 的输出是不同的，它们是矩阵。 我用 30 位小数打印它们 A和x在MATLAB和C中是一样的

A = [0, 0.100000000000000005551115123126,0,0,0;
-1705.367199390822406712686643004417 -13.764624913971095665488064696547 245874.405372532171895727515220642090 0.000000000000000000000000000000 902078.458362009725533425807952880859； 
0, 0, 0, 0.100000000000000005551115123126, 0;
2.811622989796986438193471258273, 0, -572.221510883482778808684088289738, -0.048911651728553134921284595293 ,0;
0, 0, -0.100000000000000005551115123126 0, 0]

x = [0.071662614269441649028635765717 ;
45.870073568955461951190955005586;
0.000002088948888569741376840423;
0.002299524406171214990085571728;
0.000098982102875767145086331744]

但是A*x的结果不一样，MATLAB中第二个元素是-663.792187417201375865261070430279，C中是 -663.792187417201489552098792046309

MATLAB
A*x = [ 4.587007356895546728026147320634
  -663.792187417201375865261070430279
  0.000229952440617121520692600622
  0.200180438762844026268084007825
  -0.000000208894888856974158859866];

C
A*x = [4.587007356895546728026147320634
 -663.792187417201489552098792046309
  0.000229952440617121520692600622
  0.200180438762844026268084007825
  -0.000000208894888856974158859866];

虽然差异很小，但我需要这个结果来做有限差分，此时结果会更明显

有人知道为什么吗？

Answer 1

你认为你需要多少位数？每个数字的前 16 位数字相同，这是 double 通常可以在内部表示和存储的近似数据量。您无法获得更多，即使您强制打印程序打印更多数字，它们也会打印垃圾。发生的事情是，您说过要在打印例程中使用 120 位数字……他们会打印这些数字，通常乘以余数（无论它可以是什么）由于数字以基数 2 表示，因此您通常不会得到零一旦通过数字的内部精度......一旦您的数字中没有更多位表示，打印实现不必就打印的数字达成一致。

假设您有一个只有 10 位精度的手动计算器。你会得到 120 位的数字。您开始计算并仅获得 10 位数字的结果……但您被要求打印一份包含 120 位数字结果的报告。嗯....因为整体计算不能超过 10 位，你能做什么？您使用的计算器无法为您提供所要求的位数......而且，52 位有效数字中以 10 为基数的数字不是整数位数（52 位有效数字中有 15.65355977452702215111442252567364 十进制数字） .你能做什么，你可以用零填充（很可能是不正确的）你可以用垃圾填充那些地方（这永远不会影响最终的 10 位数结果）或者你可以去 Radio Shack 并购买一个 120 位数的计算器。浮点打印例程使用计数器来指定进入循环并获得另一个数字的次数，它们通常在计数器达到其限制时停止，但不要做任何额外的努力来知道你是否发疯并指定了一个大量的数字......如果你要求 600 位，你只会得到 600 次循环迭代，但数字会是假的。

您应该期望2^52 在double 数字中存在一个差异，因为这些是用于有效数字的二进制数字的数量（这大约是2,220446049250313080847263336181641e-16，因此您必须将此数字乘以你有输出来查看舍入误差有多大，大约）如果你将你的数字乘以 663.792187417201375865261070430279，你会得到 1.473914740073748177152126604805902e-13，这是对数字中最后一个有效数字的估计。由于进行单元计算所需的大量乘法和总和，误差估计可能会大得多。无论如何，1.0e-13 的分辨率非常好（亚原子差异，值应该是长度和单位为米）。

编辑

例如，只要考虑以下程序：

#include <stdio.h>
int main()
{
        printf("%.156f\n", 0.1);
}

如果你运行它，你会得到：

0.100000000000000005551115123125782702118158340454101562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

这确实是机器可以用基数 2 浮点表示的数字 0.1 的内部表示的最（精确）近似值（0.1 恰好是一个周期数，当以基数 2 表示时）它表示是：

0.0001100110011(0011)*

所以它不能用 52 位精确地表示，无限期地重复模式 1100。在某些时候你必须削减，printf 例程继续向右添加零，直到它达到上面的表示（以 2 为底的所有有限位数都可以表示为以 10 为底的有限位数，但是converse 不成立，（因为 2 的所有因数都在 10 中，但不是 10 的所有因数都在 2 中）。

如果将0.1 和0.1000000000000000055511151231257827021181583404541015625 之间的差除以0.1，您将得到5.55111512312578270211815834045414e-17，它大约是限制1/2^52 的1/2^54 或四分之一（大约五分之一）我在上面给你看了。它是与数字0.1 最接近的可用 52 位表示的数字