〜 微分積分学演習II,複素関数論演習の数値計算レポートに 取り組んでいる人へ 〜
数値計算は便利なものだ. しかし,慣れるまでの最初のうちは, 山のようなエラーに遭遇して 嫌になることがあるかもしれない. 数値計算ができなくて困っている人は, 以下のことを試してみて欲しい. なお,以下に記す内容は 「エラーの嵐から脱出し,正しい計算結果を得るための道筋」に関する 私の経験則である. 数値計算法やfortran文法については, 数値計算の講義資料などを参照されたい.
目次苦労してプログラムを書き終わったあと,普通 まず最初に遭遇するのはコンパイルエラーだ. ここでめげてはいけない. まずは エラーメッセージが何を言っているのかちゃんと読む ことが必要. コンパイルエラーには様々な種類があるが, エラーメッセージが伝えようとしていることは,要するに
ということだ. エラーメッセージには,大概 プログラムのどこがどう間違っているのか書かれている. これを読み解けば,プログラムをどうなおせばよいかわかる訳だ.
メッセージは英語で書かれている.まずは わからない単語を全部辞書で調べる こと. 耳慣れない単語が多くて最初は戸惑うだろうし, 単語を全部調べても メッセージの内容がわからない場合だってあるだろう. でも,「どうやらこのあたりに文句を言われているようだ」 というあたりがつけばしめたもの. 怪しいと思う箇所をなおしながら試行錯誤することで, 間違いの箇所はたいていわかる.
例として, 2つの実変数 x, y に対して a=2x+y を出力するプログラムを考える.
1 program example1 2 implicit none 3 real::x,y,a 4 5 x=1.0 6 y=2.0 7 a=2x+y 8 write(*,*) 'x=',x,' y=',y 9 write(*,*) 'answer=',answer 10 stop 11 12 end program example1
このプログラムをコンパイルすると, 次のエラーメッセージが出る (コンパイラによって文言は多少違うかもしれない). 青字の部分は,エラーが出たところの行番号を示している.
0007) a=2x+y Variable X follows another operand 0009) write(*,*) ' answer=',answer ANSWER must appear in a type declaration
7行目のエラーメッセージは "変数 X が別の被演算子の次に来ている(続いている)" と言っている. もしも「被演算子」の意味がわからなくても variableとfollowの意味さえわかれば, 「変数 x が何かの後に来るのがいけないらしい (-> xの前が怪しそう)」 ということは想像がつく. 7行目の「a=2x+y」を「a=2*x+y」になおすとこのエラーはでなくなる.
9行目のエラーメッセージは "ANSWERは型宣言に現れねばいかん" と言っている. 宣言について文句を言われているようなので 9行目と宣言文(3行目)を見かえしてみると, 9行目の answer という変数が3行目で宣言されていないことがわかる. 3, 7行目の「a」を「answer」にかえるか, あるいは9行目を 「write(*,*) ' answer=',a」 にかえるかすれば,エラーは出なくなる.
エラーメッセージに指摘された行を調べると 文法ミスが見つかることが多い. ただし,上の例の9行目エラーのように, エラーの源が指摘された行よりはるか 前に存在することもある. エラーの源がわからないときは, 指摘された行より前にある文も疑ってかかる とよい.
例1ではエラーは2個しか出なかったが,長いプログラムの場合は 100個くらいのエラーがどばどば出ることもある. そういうとき,後ろの方のエラーは前のエラー (行番号が若いもの) の 副作用として出ていることが多い. したがって, エラーの嵐に出会ったときは上のエラーから順番につぶしていく と効率がよい.
コンパイルがエラーなく通れば,実行ファイル (Windowsなら *.exe) ができる. これを実行すれば計算結果が得られるはずだが, 一発でうまく実行できず,エラーが出て計算が止まることもある. これはよくあることだ. こんなことで数値計算をあきらめてはいけない. 実行時エラーが出た場合は,まずは エラーメッセージをちゃんと読む こと. もしもエラーメッセージをちゃんと読んでも なぜエラーになるのかわからない場合は, write文を使ってどこまで計算が進んでいるのか確認するとよい.
再び例で説明する.次の例2を実行すると, 実行時エラーとして「Floating point divide by zero」という メッセージが出て, 計算結果は何も出力されない. このエラーメッセージは「ゼロ割り」を指摘しているのだが, 以下ではこのメッセージの内容がわからない場合の対応を説明する.
1 program example2 2 implicit none 3 real::x,y,a,b,answer 4 5 x=1.0 6 y=2.0 7 a=2*x+y 9 b=1.0/(a**2-16.0) 10 answer=a*b 11 write(*,*) 'x=',x,' y=',y 12 write(*,*) 'Answer=',answer 13 stop 14 15 end program example2
以下の例では,エラーの原因を探すために 各行の間にwrite文を挿入した.
1 program example3 2 implicit none 3 real::x,y,a,b,answer 4 5 x=1.0 6 y=2.0 7 write(*,*) 'ここまではOKかな? (1) x=', x, ' y=',y 8 a=2*x+y 9 write(*,*) 'ここまではOKかな? (2) a=',a 10 b=1.0/(a**2-16.0) 11 write(*,*) 'ここまではOKかな? (3) b=',b 12 answer=a*b 13 write(*,*) 'x=',x,' y=',y 14 write(*,*) 'Answer=',answer 15 stop 16 17 end program example3
例3を実行すると,例2と同様に実行時エラーが出て計算は止まる. ただし,出力として次のものが得られる.
ここまではOKかな? (1) x=1.00000 y=2.00000 ここまではOKかな? (2) a=4.00000
7,9行目のwriteは出力されたが, 11行目のwriteの結果「ここまではOKかな? (3) b=...」が出力されない. このことから,9行目までは計算が進んでいて, 10行目か11行目のどちらかで計算が止まったのだということがわかる. このプログラムでは10行目で 1 を a**2-16.0=0 で割っている. これがエラーの原因である.
例3のように write文を使って次の2点を確認 すれば,大概の実行時エラーの原因はわかる.
ただし,いつもこの方法でうまくいく訳ではない. やっかいなのは,実行時エラーが出たり出なかったりする場合だ. よくあることとして 「write文がないときは実行時エラーが出るのに, write文を追加すると実行時エラーが出なくなってしまう」 現象がある. こういうときには subroutine間でのデータの受け渡しのところで 変なミスをしているとか, 配列をどこかで壊しているとか, そんな類の間違いをしていることが多い. このような間違いはwrite文の周辺を探しても発見できない. プログラム全体を疑ってかかる必要がある.
エラーの源を探す際には, コンパイラオプションを活用するのも有効だ. コンパイルするときに使う命令(f77, g77, ftn95, fortなど)には, いろんなオプションが用意されている. 「警告を沢山だしてくれ (文法ミスじゃなくても怪しそうな所は教えてくれ)」とか, コンパイラによっては「実行時に配列を壊したときは教えてくれ」とか, そういったことをお願いできる機能がある. オプションの機能と使い方はコンパイラによって異なる. 使ってみようという人は各自コンパイラの help 等で調べよ.
コンパイルエラーもなくなった,実行時エラーもなくなった... ここでよかったよかったと安心して そのままレポートを提出してはいけない. なぜなら,
から. 数値計算の間違いをなくす作業は, エラーがなくなった段階から始まると言ってもよい.
次の例4は,入力 x に対して x の3乗根 (1/3乗)を出力するプログラムである. この例ではコンパイルエラーも実行時エラーも出ない. しかし,計算結果は間違っている.
1 program example4 2 implicit none 3 real::x,answer 4 5 x=8.0 6 answer=x**(1/3) 7 write(*,*) 'x=',x,' answer=',answer 8 stop 9 10 end program example4 実行結果 x=8.00000 answer=1.00000
計算結果が間違っているという表現は正しくないかもしれない. 計算機はプログラムに書かれた命令どおり正しく計算をしている. しかしその計算は人間の意図したものとは異なっている. 人間の意図通り計算すれば,answerは 「x**(1/3)=8**(1/3)=2」となるはずだが, 計算機の出す答えは「x**(1/3)=8**(0)=1」である. この例のように, エラーはでないが プログラムには間違いがある というのはよくあることだ.
上の例は数値計算しなくても正解がわかる問題だったので, 計算結果が間違っていることはすぐにわかった. しかし,普通は数値計算で扱うのは答えがわからない問題である. 答えがわからない問題について 計算結果が正しいことを確認するのはそう簡単ではない. 難しい作業ではあるが, 計算結果が正しいかどうか確認することは プログラムを書くこと以上に大切な作業 である.
計算結果が正しいことを確認するための万能な方法はない. 問題に応じてもし可能なら,次のことを試すとよい.
より簡単な問題,答えが解析的に得られる問題などを計算してみて, 結果が正しいかどうか確認する. 例えば,曲線の長さを計算する問題なら, 直線や円弧などの計算を試す. 或いは,非線形方程式を計算する問題なら, 非線形項を0とおいて 線形方程式が正しく解けるかどうか試す.
例えば,行列 A の逆行列 B を求める問題なら, 得られた答え B をもとの行列 A にかけてみて, 単位行列が得られるかどうか確認する. 方程式 f(x)=0 を解く問題なら, 得られた答え x を f(x) に代入して 0 になるかどうか確認する.
「正しい答えならこうなるはずだ」と 自分の思うところを試してみることも必要. 例えば,力学の問題でエネルギー等の保存量が存在するならば, 保存されるべき量が一定になっているかどうか確認する.
例えば,曲線の長さを計算する問題には 曲線の分割数 N というパラメタが存在する. 常微分方程式 dx/dt=f(x,t) を数値積分する問題であれば, 時間刻み 冲 というパラメタが存在する. このようなパラメタについては,値を幾つか変えて計算し, 結果がどのように変わるか(変わらないか)確認する.
「プログラムに間違いがあるが, ある特定の入力のときだけはたまたま正しく計算できてしまう」 ということがある. 複数の入力を試すとこのような間違いが見つかることがある. 例えば,連立1次方程式 Ax=bを解く問題なら, 行列 A やベクトル b として複数のものを与えてみる. FFTのプログラムなら, 入力として sin x だけでなく,sin x + 2 cos 3 x なども与えてみる.
計算結果はどうも正しくないようだ... プログラムにはどこかに間違いがあるみたいだ... でもどこが間違っているのかわからない. そんなとき,私なら次のことを試す. 間違いを発見するための万能な方法は無い. 基本となる 出発点は全てを疑ってかかること だ. 「今目の前にあるプログラムを書いた奴 (=過去の自分) は, 自分よりも100倍まぬけでおっちょこちょいで信用できない奴なのだ, 考えられんようなお馬鹿な間違いがどこに潜んでいても不思議はない」 そう思ってプログラムを見返すと,間違いが見つかりやすいような気がする.
実行時エラーの源を見つけるときと同じように write文を使う. 計算の各段階において変数に入っている値が 自分の意図の通りになっているかどうか確認する.
自分が計算機になったつもりで,プログラムの指示にしたがって 一行一行計算結果を書き下してみる. do loopについては,i=1のとき... i=2のとき... 全て書き下してみる. if文については,条件式やそこに含まれる変数の中身を書き下す.
この作業をするときは馬鹿になることが必要. 賢い知能 (想像力) を使って 「変数 x にはこの段階では x=2 が入っているはずだから...」 などと考えてはいけない. 計算機がやるのと同じように, プログラムの命令を 上から順番に 逐一追うこと.
例えば,N×Nの行列の問題を扱っているならば, N=2とか3くらいの小さいサイズの問題を考えてみる. 微分方程式 dx/dt=f(x,t) を数値積分する問題であれば, 時間刻み 冲 を1ステップだけ進める問題を考える.
subroutineやfunctionなど各部品ごとに 単体で入出力を確認する.
「プログラムには間違いがなくて, そもそも計算させようとしていた 式自体が間違っていた」ということもある. 自分のノートやメモ書きなどが正しいかどうかを確認する.
「プログラムには間違いがなくて, 間違いがないかどうか試すテスト(検算等)の方が間違っていた」ということもある. テスト方法とその結果を改めて吟味する.
丸一日プログラムを見返しても間違いが見つからないときは, 風呂に入って一息いれてみる. 一度頭をからっぽにしてからプログラムを見返すと, 思いのほか単純な間違いが見つかることがある.