高精度数值计算

此示例说明如何通过 Symbolic Math Toolbox™ 使用可变精度算术进行高精度计算。

搜索表示接近整数的公式。例如，将 exp(sqrt(163)*pi) 计算到 30 位数。结果显示为有舍入误差的整数。

disp("Setting precision to 30 digits")
digits(30);
f = exp(sqrt(sym(163))*sym(pi));
vpa(f)

Setting precision to 30 digits
 
ans =
 
262537412640768743.999999999999
 

当您将相同的值计算到 40 位数时，可以看到该值实际上不是整数。

disp("Setting precision to 40 digits")
digits(40);
vpa(f)

Setting precision to 40 digits
 
ans =
 
262537412640768743.9999999999992500725972
 

接下来，假设有一系列数字，其值最高达到 exp(1000)。确定小数点后的正确位数，以便正确计算这些数字。找到上界 exp(1000) 所需的最小工作精度。

disp("Compute the required working precision")
disp(">> d = log10(exp(vpa(1000)))")
d = log10(exp(vpa(1000)))

Compute the required working precision
>> d = log10(exp(vpa(1000)))
 
d =
 
434.2944819032518276511289189166050822944
 

请先设置所需的精度，然后再开始调用依赖该精度的函数。例如，使用函数 round、vpa 或 double。

digits(ceil(d) + 50);

现在，假设有一系列 exp(sqrt(n)*pi) 形式的数字，n 的值为从 1 到 1000。检查此形式的数字是否接近某个整数。您可以从它们的小数部分的直方图中看到这一点。

A = exp(pi*sqrt(vpa(1:1000)));
B = A-round(A);
histogram(double(B), 50)

计算是否有 exp(n) 形式的近整数。

A = exp(vpa(1:1000));
B = A-round(A);
find(abs(B) < 1/1000);

现在您可以看到 A 的元素的小数部分分布相当均匀。

histogram(double(B), 50)