可能な限りシンボリック計算を使用

計算をシンボリックに行うことが推奨されるのは、厳密にシンボリック計算には丸め誤差が発生しないためです。たとえば、以下のような標準的な数学定数には、Symbolic Math Toolbox™ でのシンボリックな表現が存在します。

pi
sym(pi)

ans =
    3.1416

ans =
pi

不必要な数値近似の使用は避けるようにします。浮動小数点の数字は定数を近似しますが、それ自体は定数ではありません。このような近似を使用すると、結果が不正確になる可能性があります。たとえば、特殊関数 heaviside は、sym(pi) の正弦と pi の数値近似の正弦とで以下のように異なる結果を返します。

heaviside(sin(sym(pi)))
heaviside(sin(pi))

ans =
1/2

ans =
     1

精度を上げて計算を実行

リーマン予想では、ゼータ関数 ζ(z) のすべての非自明な零点は同一の実数部 ℜ(z) = 1/2 を有していると提唱されています。ゼータ関数の零点を見つけるには、その絶対値 |ζ(1/2 + iy)| をプロットします。次のプロットは、ゼータ関数 |ζ(1/2 + iy)| の最初の 3 つの非自明な根を示しています。

syms y
fplot(abs(zeta(1/2 + i*y)), [0 30])

数値ソルバー vpasolve を使って、このゼータ関数に存在する最初の 3 つの零点を近似します。

vpasolve(zeta(1/2 + i*y), y, 15)
vpasolve(zeta(1/2 + i*y), y, 20)
vpasolve(zeta(1/2 + i*y), y, 25)

ans =
14.134725141734693790457251983562
 
ans =
21.022039638771554992628479593897
 
ans =
25.010857580145688763213790992563

次に、同じ関数の実数部を少し大きくし、$$\zeta \left(\frac{1000000001}{2000000000}+iy\right)$$ とします。リーマン予想によると、この関数ではどの実数値 y に対しても零点がありません。有効小数桁数 10 桁で vpasolve を使用すると、ソルバーは以下のようにゼータ関数の (存在しない) 零点を見つけます。

old = digits;
digits(10)
vpasolve(zeta(1000000001/2000000000 + i*y), y, 15)

ans =
14.13472514

桁数を増やすと、この結果が不正確であることがわかります。ゼータ関数 $$\zeta \left(\frac{1000000001}{2000000000}+iy\right)$$ には、実数値 14 < y < 15 に零点が存在しません。

digits(15)
vpasolve(zeta(1000000001/2000000000 + i*y), y, 15)
digits(old)

ans =
14.1347251417347 + 0.000000000499989207306345i

計算を続けるため既定の桁数に戻します。

digits(old)

シンボリックと数値のそれぞれの結果を比較

インデックスが半整数のベッセル関数は、厳密なシンボリック式を返します。このような式を浮動小数点数で近似すると、結果が不安定になることがあります。たとえば、次のベッセル関数の厳密なシンボリック式は以下のとおりです。

B = besselj(53/2, sym(pi))

B =
(351*2^(1/2)*(119409675/pi^4 - 20300/pi^2 - 315241542000/pi^6...
 + 445475704038750/pi^8 - 366812794263762000/pi^10 +...
 182947881139051297500/pi^12 - 55720697512636766610000/pi^14...
 + 10174148683695239020903125/pi^16 - 1060253389142977540073062500/pi^18...
 + 57306695683177936040949028125/pi^20 - 1331871030107060331702688875000/pi^22...
 + 8490677816932509614604641578125/pi^24 + 1))/pi^2

この式を、vpa を使用して 10 桁の精度で近似します。

vpa(B, 10)

ans =
-2854.225191

次に、浮動小数点のパラメーターを指定してベッセル関数を呼び出します。この 2 つの近似の有意差は、どちらか一方もしくは両方の結果が正しくないことを示しています。

besselj(53/2, pi)

ans =
   6.9001e-23

数値処理の精度を上げ、より正確に B を近似します。

vpa(B, 50)

ans =
0.000000000000000000000069001456069172842068862232841396473796597233761161

関数または式をプロット

結果をプロットすると、不正確な近似を認識するのに役立ちます。たとえば、次のようにベッセル関数を数値近似すると、以下が返されます。

B = besselj(53/2, sym(pi));
vpa(B, 10)

ans =
-2854.225191

このベッセル関数を 53/2 周辺の x 値についてプロットします。この関数のプロットを見ると、近似が不正確であることがわかります。

syms x
fplot(besselj(x, sym(pi)), [26 27])