MATLAB 使用 vpa() 时 mod() 的奇怪行为答案

【问题标题】：MATLAB strange behaviour with mod() when using vpa()MATLAB 使用 vpa() 时 mod() 的奇怪行为
【发布时间】：2020-07-08 23:58:10
【问题描述】：

如果我们有以下代码：

base = vpa(1); 
height = vpa(2.2);
mod(2*base + height + height, 2 * (base + height))

这会产生 6.4 的输出。我希望结果为 0，而数值解确实为 0。但我需要将符号值与 vpa() 一起使用。

我做了一些实验来找出原因并发现：

simplify(2*base + height + height < 6.4)
simplify(2 * (base + height) == 6.4)

两者都给TRUE。所以相同的（数字）表达式更小，等于 6.4。

我应该怎么做才能解决这个问题并得到 0 的答案？是什么导致了这个问题？

【问题讨论】：

这确实是一个意外的行为。解析上相同的表达式在符号数学工具箱中并不相同。我猜，它与舍入有关，但我不知道符号计算是如何工作的。

标签： matlab symbolic-math

【解决方案1】：

问题在于vpa 提供了任意大的精度，但并不精确。首先，注意vpa 与第二个输入使用digits 精度数字，默认为32。当你这样做时

height = vpa(2.2)

这是一样的

height = vpa(2.2, 32)

假设digits 是32。所以height 会有32 的精度数字，但不会是精确的。要看到这一点，请更精确地评估定义的height：

>> vpa(height, 32)
ans =
2.2
>> vpa(height, 40)
ans =
2.2
>> vpa(height, 50)
ans =
2.2000000000000000000000000000000000000001469367939

这种不准确性导致2*base + height + height 和2 * (base + height) 之间存在数值差异，两者实际上都不等于6.4：

>> base = vpa(1); 
>> height = vpa(2.2);
>> vpa(2*base + height + height, 50)
ans =
6.3999999999999999999999999999999999999988245056492
>> vpa(2 * (base + height), 50)
ans =
6.4000000000000000000000000000000000000002938735877

因此，即使mod(2*base + height + height, 2 * (base + height)) 似乎是0

>> mod(2*base + height + height, 2 * (base + height))
ans =
6.4

but it's _not_:

>> vpa(mod(2*base + height + height, 2 * (base + height)), 50)
ans =
6.3999999999999999999999999999999999999988245056492

请注意，后一个结果中与6.4 的偏差不等于上述两个小偏差的总和；相反，它等于第一个。数值误差不保证是累加的。

简而言之，vpa 减少但不能完全避免数值精度误差。

如果我们增加vpa 的位数会怎样？这会提供更高的精度并可能解决问题吗？

>> base = vpa(1, 1000); 
>> height = vpa(2.2, 1000);

奇怪的是，尽管使用了更多的数字，但我们还是得到了和以前一样的不准确：

>> vpa(2*base + height + height, 50)
ans =
6.3999999999999999999999999999999999999988245056492
>> vpa(2 * (base + height), 50)
ans =
6.4000000000000000000000000000000000000002938735877
>> vpa(mod(2*base + height + height, 2 * (base + height)), 50)
ans =
6.3999999999999999999999999999999999999988245056492

因此，避免精度错误的唯一方法是将vpa 替换为符号变量，它们是精确：

>> base = sym(1)
base =
1
>> height = sym(2.2)
height =
11/5

关于如何定义符号变量的题外话是有序的。请注意，sym(2.2) 有一个潜在问题，它首先将2.2 定义为double 浮点数，具有其固有的不准确性，然后然后被转换到sym。在这种情况下，这不是问题，因为2.2 的浮点数字表示恰好是精确的。事实上，我们可以检查 Matlab 是否显示 height = sym(2.2)，这是准确的。此外，即使 double 表示不准确，Matlab 也会尝试猜测您的意图，并且通常会成功：

>> sym(3.141592653589793)
ans =
pi

（Matlab 假设我们指的是数字 pi）...但并非总是如此：

>> sym(sqrt(777))
ans =
777^(1/2)
>> sym(sqrt(777777))
ans =
7757421003204227/8796093022208

(sqrt(777) 给出double 结果27.874719729532707，它被sym 识别为777 的平方根的近似值。另一方面，sqrt(777777) 给出8.819166627295348e+02，它未被sym 识别为777777 的平方根的近似值。

>> sym(7777)
ans =
7777
>> sym(7777777777777777777)
ans =
7777777777777777664

（7777 完全表示为double，但7777777777777777777 不是，因为它超过了2^53。）

因此，可以肯定的是，定义符号变量的安全方法是只使用小整数、其他符号变量或精确的字符串表示：

>> height = sym('22/10')
height =
11/5

现在，将base 和height 正确定义为符号变量，我们得到了预期的结果：

>> mod(2*base + height + height, 2 * (base + height))
ans =
0
>> vpa(mod(2*base + height + height, 2 * (base + height)), 1000)
ans =
0.0

【讨论】：

有趣！我想你可以做height = vpa(22, 1000) / vpa(10, 1000)？两者都应该准确表示？
@CrisLuengo 这有同样的问题：结果是2.2000000000000000000000000000000000000001469...。不知何故，2.2 和 vpa 是不准确的，即使它与 double 完全一样

【解决方案2】：

mod(2*base + height + height, base+base+height+height)

这行代码生成 0。仍然不确定究竟是什么导致了幕后问题。 product 和 sum 之间似乎有明显的区别，因为下面的代码行是 FALSE

mod(base+base+height+height, 2*(base+height))

【讨论】：