「二次方程式の解を求めるプログラムのレビュー」で,実数解の場合に解の公式だけで解を求めてはならないと述べたが,Python の sympy がどうやっているのかレビューした結果を示す。
以下の二次方程式を解け。
x**2 +1.0000000000000001*x + 0.0000000000000001 = 0
結論から述べると,正しいやり方は「二次方程式の係数を Rational() で与え,solve() で解を求める」ということである。
x**2 +1.0000000000000001*x + 0.0000000000000001 = 0
結論から述べると,正しいやり方は「二次方程式の係数を Rational() で与え,solve() で解を求める」ということである。
>>> from sympy import *
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)
係数を Rational() で定義する。
>>> var('a:z')
(a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)
係数を Rational() で定義する。
>>> b = Rational(10000000000000001, 10000000000000000)
>>> c = Rational( 1, 10000000000000000)
>>> print(b, c)
10000000000000001/10000000000000000 1/10000000000000000
Rational() ならば精度は完全に保証される。
>>> print(b.evalf(30))
1.00000000000000010000000000000
>>> print(c.evalf(30))
1.00000000000000000000000000000e-16
solve() で解を求める
>>> ans = solve(x**2 + b * x + c)
>>> print(ans)
[-1, -1/10000000000000000]
つまり,
>>> print(ans[0].evalf(30))
-1.00000000000000000000000000000
>>> print(ans[1].evalf(30))
-1.00000000000000000000000000000e-16
>>> print(ans[0].evalf(30))
-1.00000000000000000000000000000
>>> print(ans[1].evalf(30))
-1.00000000000000000000000000000e-16
以下は読まなくて良い。
ダメなやり方
python の float を Float() で変換しても,evalf() で有効数字を 30 桁表示させれば精度が足りないのは明らかである。
>>> b = Float(1.0000000000000001)
>>> c = Float(0.0000000000000001)
>>> print(b.evalf(30), c.evalf(30))
1.00000000000000000000000000000 9.99999999999999979097786724035e-17
文字型で指定してやれば,精度は増す。
>>> b = Float('1.00000000000000010000000000000000000000000000')
>>> c = Float('0.00000000000000010000000000000000000000000000')
>>> print(b.evalf(30), c.evalf(30))
1.00000000000000010000000000000 1.00000000000000000000000000000e-16
十分(?)な精度の b, c で二次方程式を解くと以下のようになる。
>>> ans = solve(x**2 + b * x + c, x)
デフォルトでの表示は正確な解が得られたと見まがう。
>>> print(ans)
[-1.00000000000000, -1.00000000000000e-16]
しかし,個別の解を .evalf(30) で表示させると,精度が足りないことが明白になる。
[-1.00000000000000, -1.00000000000000e-16]
しかし,個別の解を .evalf(30) で表示させると,精度が足りないことが明白になる。
>>> print(ans[0].evalf(30))
-0.999999999999999888977697537484
>>> print(ans[1].evalf(30))
-1.00000000000000010235730316482e-16
解の公式だけで実数解を求める場合には,大小どちらの解も精度は十分であることがわかる。(何故 solve() の場合と違うの か?)
>>> x1 = (-b - sqrt(b**2 - 4 * c)) / 2
>>> print(x1.evalf(30))
-1.00000000000000000000000000000
>>> x2 = (-b + sqrt(b**2 - 4 * c)) / 2
>>> print(x2.evalf(30))
-1.00000000000000000000000000000e-16
>>> print(x1.evalf(30))
-1.00000000000000000000000000000
>>> x2 = (-b + sqrt(b**2 - 4 * c)) / 2
>>> print(x2.evalf(30))
-1.00000000000000000000000000000e-16
※コメント投稿者のブログIDはブログ作成者のみに通知されます