実習課題1: 指数関数型の時間発展

変数 S の単位時間あたりの変化がS自体に比例する問題を考える。 比例定数を a とする。 (6月2日の講義では減衰問題に限ったため比例定数を -a (aが正)としたが、 成長問題も含めるにはこのほうがわかりやすいので変更した。 タンクモデルの話と対応づける場合には符号に注意してほしい。)

dS/dt = a S ... (1)

S(t) = S(0) exp(a t) ... (2)

これを、前方差分 [差分法の教材参照]で計算してみる。 比較のため、解析解を計算して同様な形で書き出すプログラムも用意した。

ソースプログラム (表の中からリンクされています)

前方差分

解析解 (言語コンパイラ・インタプリタ常備のexp関数を使う)

プログラムへの値の入力

次の数値を、 スペース区切り(連続するスペースは1つのスペースと同じとみなす)で 同じ行に順に与える。

szero (浮動小数点数)

S(0)、時刻t = 0でのSの値。

a (浮動小数点数)

式(1)のaの値。

deltat (浮動小数点数)

Δt、差分の時間刻み。

nstep (整数)

計算するステップ数。計算の終わりの(モデルの世界での)時刻は nstep*deltat となる。

intvp (整数)

計算するステップのうちいくつおきに結果を書き出すか。 毎回書き出す場合は1とする。

注

• 実用的には、nstep, intvpを指定するよりも、 計算の終わりの時刻、結果を書き出す時間間隔を、 それぞれモデルの世界の時間の単位で指定するようにしたほうがわかりやすい。 しかし、時間間隔からステップ数を求める計算を単純にプログラミングすると、 浮動小数点数どうしの割り算の結果は本来整数になるべきものでも ぴったり整数になるとは限らないので、 切り捨てによって意図したものと違う値になってしまう心配がある。 これの対策をするとプログラムが長くなるのでここでは避けた。
• 解析解の場合、nstep*deltat, intvp*deltatのそれぞれの積には意味があるが、 nstep, intvp, deltatというふうに分ける意味がなくなっている。 しかし便宜上、差分の場合と同じ形にした。
• Awk言語インタプリタには、データを与えたあと、 データの終了を知らせてやる必要がある。 これにはコントロールキーを押しながらDを押せばよい(はずである)。
• C言語の例文では、入力変数を読むのに関数scanfを使った。 これはちょっと危険だが [別記事参照]、 プログラムを短くするために使った。
• 数量の単位は数値計算向きプログラム言語では直接扱わないので、 ユーザーが意識しなければならない。 今回の例文では、時間の単位、Sの単位はいずれもユーザーが自由に決められる。 ただし a の単位は自動的に「1/(時間の単位)」になる。 一般には、 物理法則の式は、時間の単位が秒であることを前提としている場合が多いので、 違う時間の単位を使う場合は、秒との間の換算をする必要があることがある。

出力

結果は標準出力に書くようにした。 とくに指定しないと画面に表示されるが、その形ではあまり実用的でない。 リダイレクション(redirection)機能 [別記事参照] を使ってファイルに書き出すことを想定している。

結果の主要部分は、 各行に、時刻tと、その時刻のSの値の、2つの数値が書かれていて、 そのような行が、書き出し対象の時間ステップの数だけならんでいるものである。 スペース区切り2列のデータとも言える。

主要部分の前に、プログラム名と、入力変数の名前・値を書きだしておく。 これらの行は、初めに「#」の字を入れておく。 「#」で始まる行はgnuplotではコメント扱いされるので、 この行が含まれたままのファイルを作図することができる。 プログラム名を書き出しているところの文字列は仮のものなので、 各自の判断で変更してほしい。 入力変数の名前・値を書きだすところでは数量の単位も書くようにするとよい。

使い方の例

Awk

awk -f ex1.awk > ex1out1.txt
100  -0.1   0.1   100  10
control-D 

C

gcc ex1.c
./a.out > ex1out1.txt
100  -0.1   0.1   100  10

Fortran

g77 ex1.f
./a.out > ex1out1.txt
100  -0.1   0.1   100  10

注

• CとFortranの場合の a.outはLinuxなどUnix系OSの場合の実行ファイル名。 CygWinの場合は、a.exeとなるので、 それを実行するには上記のa.outの代わりにaで よいと思われる(ただし未確認)。
• Cで、expfunc.cをコンパイルする際には ライブラリ指定「-lm」で 関数expを含むmath libraryをリンクする必要がある。 今後の演習問題では、 差分の場合にもこれが必要になる場合もあるかもしれない。
• ここで使ったFortran文法は、 基本的にはFortran 77 (1977年版標準)だが、 Fortran 90で採用されたいくつかの構文 (DO...END DO, DO WHILE...END DO, DOループから脱出するEXIT文) を加えたものである。

作図

計算の結果をグラフにしてみることを勧める。 グラフ作成ソフトウェアとしては、gnuplotを勧める。 上に述べた実行例のようにファイルを作ったとすれば、 gnuplotの中で

plot "ex1out1.txt"

のようにすれば、横軸に時刻、縦軸にSの値をとったグラフができる。 plotコマンドには、ファイル名をコンマでくぎって複数ならべる ことができ、そうすると複数のデータのグラフが同じ画面に表示される。 [gnuplot教材ページ] も参考にしてほしい。