2012年1月16日月曜日

数値計算に勝つ!

いつからだろう.Excelや関数電卓を使う機会が増えたのは.

数値計算はとても便利!
だけど,よくわからないで使っていると痛い目にあう.
答えが間違っていたり,答えがでなかったり.

最近そんなことが起こり始めたので,その失敗談(?)と原因を考察してみた.

■Case1
4次ルンゲクッタ法による微分方程式の数値計算.

工学系の数理という授業の課題で
「粘性,圧縮性および表面張力を無視できる完全流体中の単一球形気泡の振動」
の微分方程式


(この式の導出考えるのおもしろそうだなぁ.今は時間がねえ.)
を線形近似した解とルンゲクッタ法による解を比べよ.

という課題がでて,tの刻みを0.00002としてルンゲクッタ法を用いると




減衰振動するんだけど,ほかの人が刻みを0.000002にしたところ減衰しないらしい.

そんなのはおかしい!
だって刻みを小さくすればするほど精度があがるのだから!!
減衰しない?
減衰項落ちてる!?
これたぶん有効数字で減衰項カットされてるじゃ!?
加工学の授業中に友達のスマートフォンでgoogleとExcel2007は15乗までしか計算しないらしい.
おう!ビンゴ!減衰項が切り捨てられている.

減衰項落ちた原因は,有効数字!


■Case2
これはhotなネタで今日,直面した


なる方程式を関数電卓を使って解くという.
これを素直に打ち込んでソルブ機能を用いると
約5秒くらいの時間経過ののち
X=0.99258・・・
という計算結果がでてくる.
今回は幸運にも答えがわかっていたため,これがおかしいということに気づく.
とりあえず
対数をとって

もしくは両辺1-Xで割って

という式に変えれば正しい答え
X=-5.8339・・・
がでる.

まぁでもこうなった原因を考えたいわけで(周りのみんなはなぜ興味ない!)
いろいろ考えた結果.
おそらく!(いまテスト前だから正しいか調べるの面倒.ブログかくなって?かきたいのっ!)
あれだ.たぶんテイラー展開とかなんらかの級数展開してるはずだから,Xe^Xとかいうやつの級数の展開は収束が遅かったりなんかして数字が切り捨てられすぎて答えがあわないんだ!
定性的にはよさそうじゃないか!
だって!
正しい答えのでる2式は級数展開簡単そうじゃん!!
答え知っている人は教えてくださいー.





数値計算はおもしろい.しっかり勉強しておくと役に立ちそうだなあ.
というかしっかり理解しないで,使われたら怖いなぁ.

2 件のコメント:

  1. 何か前提条件があるのかもわかりませんが,方程式だけみれば,X=0.99258...も誤った解ではありませんよ.f(x)=(1-x)exp(x)-0.02のグラフを描いてみてください.
    関数電卓が中で何をしてるかはわかりませんが,ふつう非線型方程式を数値的に解くアルゴリズムの代表格といえばNewton法でしょう.この問題をNewton法で解くと,初期値が0だと計算が破綻し,正だとX=0.99258...,負だと-5.8339...に収束します.これもグラフの様子とNewton法の性質からすぐわかるものです.
    変形した2式で解がX=0.99258...と求まらないのは,(1-x)が0に近いために1/(1-x)やlog(1-x)の数値計算が難しく,誤差が入った為だと思います.

    返信削除
    返信
    1. コメントどうもありがとうございます.
      コメントもらえるのはうれしいです.

      たしかにX=0.99258..も正しいですね!調べるべきでした.そうなるとぼくの予想は間違っていますね!
      Newton法だとたしかに収束する解のみ表示されることになりそうです.
      これを参考にテスト明けに詳しく調べてみることにします.

      どうもありがとうございました!

      削除