高精度の数値計算

この例では次を使用します。

この例では、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);

ここで、1 から 1000 までの n について、式 exp(sqrt(n)*pi) の数値を考えます。この式の数値が何らかの整数に近いかどうかを確認します。それらの小数部のヒストグラムプロットからこれを確認できます。

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)