Mapleでの検算の限界について

質問

Win2kでMaple7を使用して化学計算をしています。

xz/y=1.75×10^(-5)
y+z=0.05
x=z+w
wx=10^(-14)

という連立方程式でxの値を求めようとして
w,y,zを消去して
x^3+1.75^(-5)*x^2-(1.0*10^(-14)+0.05*1.75^(-5))*x-1.75^(-5)*1.0*10^(-14)=0
として
> solve({x^3+1.75^(-5)*x^2-(1.0*10^(-14)+0.05*1.75^(-5))*x-1.75^(-5)*1.0*10^(-14)=0},{x});
{x = -.9350613693e-1}, {x = -.2000000000e-12}, {x = .3257914223e-1}

となりました。
x = .3257914223e-1
という解の検算をしてみました。

> p:=x^3+1.75^(-5)*x^2-(1.0*10^(-14)+0.05*1.75^(-5))*x-1.75^(-5)*1.0*10^(-14);
　　　p := x^3+.6092699470e-1*x^2-.3046349735e-2*x-.6092699470e-15
> subs(x=0.03257914223,p);
　　　-.4060926995e-13

となりました。
これは0に近いと見なさざるえないのですかね。

回答No.1

yaksa · ベストアンサー

数値で計算すれば、誤差がでます。
普通は、数値は16桁程度の有効桁です。

せっかくMapleなんだから、記号計算すれば誤差はなくなると思います。
1.75を175/100
0.05を5/100
と書けば、記号のまま計算してくれるでしょう。
（めちゃめちゃ複雑な答えになると思いますが。）

それから、ついでに数値で計算するときには、
p := x^3+.6092699470e-1*x^2-.3046349735e-2*x-.6092699470e-15
ではなくて、
p := ((x+.6092699470e-1)*x-.3046349735e-2)*x-.6092699470e-15
の形にして計算すると、桁落ちの影響がすくなくなって、もう少し0に近づくかもしれません。
Mapleがデフォルトで整式は、こうやって計算してる可能性が高いのでおそらく変わらないとは思いますが。

補足

ご回答有難うございます。

> 数値で計算すれば、誤差がでます。
> 普通は、数値は16桁程度の有効桁です。
有難うございます。覚えておきます。

> せっかくMapleなんだから、記号計算すれば誤差はなくなると思います。
> 1.75を175/100
> 0.05を5/100
> と書けば、記号のまま計算してくれるでしょう。
> （めちゃめちゃ複雑な答えになると思いますが。）
> それから、ついでに数値で計算するときには、
> p := x^3+.6092699470e-1*x^2-.3046349735e-2*x-.6092699470e-15
> ではなくて、
> p := ((x+.6092699470e-1)*x-.3046349735e-2)*x-.6092699470e-15
> の形にして計算すると、桁落ちの影響がすくなくなって、もう少し0に近づく
> かもしれません。
> Mapleがデフォルトで整式は、こうやって計算してる可能性が高いのでおそら
> く変わらないとは思いますが。
有難うございます。
早速、試してみましたら何故か虚数が含まれる解ばかりが出力されてしまいました(下記)。
plotコマンドでグラフを出力すると(下記のように)x軸と3箇所で交わっているので虚数解は出ない筈だと思うのですが、、、
これはどうしてなんでしょうか?

後、
解を少数で出力させるにはsolveコマンドにどう記述したらいいでしょうか?

> solve({((x^2+(175/10000000)*x)-(1/100000000000000+(5/100)*(175/10000000)))*x-(175/10000000)*(1/100000000000000)=0},{x});

65632657
--------
(1/3) 15000000
{x = 1/15000000 %1 + -------- - 7/1200000}, {x =
1/3
%1

(1/3) 65632657 1
- 1/30000000 %1 - -------- ------- - 7/1200000
30000000 (1/3)
%1

/ (1/3) 131265314 1 \
+ 1/800000 I sqrt(3) |2/75 %1 - --------- -------|}, {x
| 75 (1/3)|
\ %1 /

(1/3) 65632657 1
= - 1/30000000 %1 - -------- ------- - 7/1200000
30000000 (1/3)
%1

/ (1/3) 131265314 1 \
- 1/800000 I sqrt(3) |2/75 %1 - --------- -------|}
| 75 (1/3)|
\ %1 /

%1 := -8613950975 + 12 I sqrt(1962833543397105872547)

> plot(((x^2+(175/10000000)*x)-(1/100000000000000+(5/100)*(175/10000000)))*x-(175/10000000)*(1/100000000000000),x=-0.002..0.002,y=-0.0000000005..0.0000000005);