前に述べた減衰( a < 0 )または増大( a > 0 )現象を表わす微分方程式
d Z / d t = a Zの数値解を差分で求め、グラフにしてみよう。
数値解の解析解(真の解)からの誤差の符号は、
また、誤差、つまり「数値解 - 解析解」のグラフを作ってみよう。
プログラムは、 forward.p, forward.f, forward.c である。 内容の説明は Fortran、 Pascal、 C、それぞれ別の記事に書く。 まず、例文のプログラムを自分のところにコピーして動かしてみよう。 また、Fortran の場合を例にあげるが、PascalやCでも同様である。)
まずforward.fを作業ディレクトリに持ってくる (ファイル名は変えてもよいが、Fortranの場合 .f で終わるようにする)。 % f77 forward.f % a.outすると、端末に数値が表示されるはずである。 (このうち、「#」で始まる行は、数値計算の条件を設定するために与えた定数 [プログラム言語の立場では変数であるが、数値実験の立場ではパラメータ (parameter)という]の値を表示しているものである。 計算結果、数学的に求めた真の解[の、計算機の数値の精度の範囲での近似] である。)
ファイルに書いてあとで見るためには、
% a.out > forward.out % less forward.outなどとすればよい。
(出力を書き出すファイルは単なるテキストファイルなので、 ファイル名はなんでもよいが、最後の部分がf, c, p, o, s, a, e, r, ps, tex などは、何かのプログラムで予約された意味があるので避けよう。 outはoutputの略で、forward.outは「forward.fの出したもの」のつもり。 もちろん、「a.out」という組み合わせだけは予約された意味があるので注意。)
出力結果をファイルにリダイレクト(redirect)して gnuplotでグラフ化してみよう。 gnuplotは、単に
% gnuplotで起動する。 「gnuplot>」および、継続行の場合の「>」は gnuplotが出してくる プロンプトで、キーボードから入れる文字ではない。 終了は(「gnuplot>」プロンプトに対し) quit または exit である。
たとえば、次のようにすると、差分計算の結果を点のような記号で、真の解 (時間についてはサンプリングされたものだが)を線で表示するので、 両者の系統的なずれを見ることができる。「1:2」はデータの1列めを横軸、 2列めを縦軸にとって表示することを示している。
gnuplot> plot "forward.out" using 1:2 with points, \ > "forward.out" using 1:3 with lines gnuplot> quit
(注: 上の plot 文は、コンマのあと続けて2回めの "forward.out"を 書いてしまってもよい。その場合、「\」は不要。 「\」はバックスラッシュで、キーボードによっては「￥」になっている。 これは、gnuplotで一つのコマンドを複数行にまたがって書くときに、続きが あることを示している。 また、続きの行の頭の「>」はgnuplotが出してくるもので、 キーボードから入れる文字ではない。)
ここでは2つのデータが同じファイルにあったが、複数のファイルからのデータを いっしょにプロットする場合も、
plot "ファイル名1", "ファイル名2", "ファイル名3" ...のように、ダブルクォートでくくったファイル名を列挙すればよい。 同時に使えるファイル数の限度は、実際上、gnuplot内のコマンド行の長さに 制約され、経験によれば20個程度である。
その他、gnuplotの使い方に関しては 別記事参照。
プログラムの中で、次の定数を、代入文の右辺に書いてある。 [プログラム言語の立場では左辺が変数、右辺が定数であり、 FortranのPARAMETER文は使っていないが、 数値実験の立場で、このようなものを実験のパラメータ(parameter)という]。
tzero = 0.0d0 時間発展の開始時刻 deltat = 0.1d0 時間ステップ (解析解では結果の出力間隔) nt = 100 終了までの時間ステップ数 zzero = 1.0d0 初期値 a = -1.0d0 方程式の係数 aこの値をいろいろと変えて比較してみよう。
% cp forward.f forw01.f % mule forw01.f 編集開始(エディタの例としてmuleをあげている) c-x c-s 保存 (muleのコマンド) c-x c-c muleの終了 % f77 forw01.f % a.out > forw01.out 最後にコンパイルされたものが実行されるできれば、具体的(自然・社会科学的)問題と、量の次元や単位に注意して、 関連づけてみよう。
パラメータの値を、プログラム中に定数として書くと、 変更するたびにコンパイルしなおさなければならない。 ほんとうは、これらの値は、READ文で読むほうが望ましい。 しかし、READ文の前に、何を入れたらいいのか人に伝える(プロンプトを出す) WRITE文もほしくなる。ところが、画面に出すべきプロンプトと、ファイルに redirectすることが多い計算結果の出力の行き先を区別するためには、 標準入出力以外のファイルの読み書き あるいは 標準エラー出力 を使わなくてはならない。 ここではプログラムで使う入出力機能を標準入出力だけに限ることを優先し、 妥協として、定数はプログラム中に書く方針をとることにした。
WRITE文の変数並びの変数の数はいくつでもよい。 ただし、書き出す変数並びと書式が対応している必要がある。 Fortranの場合で言えば、
WRITE(*,710) t, z, ztrue
710 FORMAT(' ',3F20.10)
を
WRITE(*,710) t, z, ztrue, z-ztrue
710 FORMAT(' ',4F20.10)
のように、write文とformat文を対応づけて変える。
上の例は、第4列に「誤差」の値を書き出すようにしたもので、
% a.out > forward.out2と、ファイルに計算結果を書き出しておき、 gnuplotで
gnuplot> plot "forward.out2" using 1:4 with pointsなどという形で作図することができる。