c - C和MATLAB中的矩阵乘法，结果不尽相同

Question

我正在使用4Rungekutta在MATLAB和C中解决DGL（Ax + Bu = x_dot），
 A是5x5，x是5x1，B 5x1，u 1x1，u是正弦函数的输出（2500点），
在第45次迭代之前，MATLAB和C中4Rungekutta的输出都是相同的，但在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的一部分有所不同，因为这些是用于有效位数的二进制数字的数量（这是aprox 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的内部表示的最大（最精确）近似值（当以基数2表示时， 0.1恰好是一个周期数）。其表示为：

0.0001100110011(0011)*

因此无法用52位无限重复地重复模式 1100来准确表示它。在某些时候，您必须剪切，然后 printf例程继续向右添加零，直到到达上面的表示形式为止（以2为底的所有有限数字都可以表示为以10为底的有限数字，但是相反不是正确的（因为2的所有因子都在10中，但并非10的所有因子都在2中）。

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