时间戳

首页 > 解决方案 > C和MATLAB中的矩阵乘法，不同的结果

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

问题描述

我正在使用 4Rungekutta 在 MATLAB 和 C 中求解 DGL(A x + B u = 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];

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

有谁知道为什么？

标签： cmatlab

解决方案

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

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

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

编辑

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

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

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

0.100000000000000005551115123125782702118158340454101562500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

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

0.0001100110011(0011)*

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

如果您将两者之间的差异分开0.1，您将得到大约 1/2^54 或我在上面向您展示的限制的大约四分之一（大约五分之一） 。它是用 52 位表示的最接近数字的数字0.10000000000000000555111512312578270211815834045410156250.15.55111512312578270211815834045414e-171/2^520.1

推荐阅读