4. f(x) の求積問題
4.1. 設計
最後も課題を確認しておきます。課題は「求積問題で適当な関数
乱数の発生については、Buffonの針と同様に、先程作成したrandコマンドがそのまま使えます。あとは乱数を読み込んで定積分するプログラムを作ることになります。今回は入門的モンテカルロ法で定積分するプログラムをcrudeコマンド、あたりはずれモンテカルロ法で定積分するプログラムをhitormissコマンドと呼ぶことにします。テキストにある通り、入門的モンテカルロ法(crudeコマンド)では一つ乱数を発生させることが1ステップに相当します。また、あたりはずれモンテカルロ法(hitormissコマンド)では二つの乱数で二次元座標上の点を表現するので、二つ乱数を発生させることが1ステップに相当します。
4.1.1. crudeコマンド
crudeコマンドで乱数以外に必要になる人為的なパラメータは、被積分関数
foo[bar]% ./rand 10 | ./crude a の値 b の値 Enter
例えば、それぞれ
foo[bar]% ./rand 10 | ./crude 0 1 Enter
.
.
.
foo[bar]%
上の例では、randコマンドが10個乱数を発生させていますので、crudeコマンドは10ステップ実行されることになります。
4.1.2. hitormissコマンド
hitormissコマンドで乱数以外に必要になる人為的なパラメータは、被積分関数
foo[bar]% ./rand 10 | ./hitormiss a の値 b の値 c の値 Enter
例えば、それぞれ
foo[bar]% ./rand 10 | ./hitormiss 0 1 1 Enter
.
.
.
foo[bar]%
上の例では、randコマンドが10個乱数を発生させていますが、hitormissコマンドは1ステップで2つの乱数をとるので、5ステップだけ実行されることになります。
4.2. 作業用ディレクトリの作成
次のようにpwdコマンドを入力して、現在自分がホームディレクトリにいることを確認してください。
foo[bar]% pwd Enter
/home/appi/foo
foo[bar]%
もし違っていたらcdコマンドでホームディレクトリに移動します。
foo[bar]% cd Enter
foo[bar]% pwd Enter(確認のため)
/home/appi/foo
foo[bar]%
求積問題のプログラムを作成するための作業用ディレクトリを次のように作ります。
foo[bar]% mkdir integral Enter
次にこのディレクトリに移ります(繰り返しになりますが、このときもtcshのファイル名補完機能を用いればタイプ量を減らすことができます)。
foo[bar]% cd integral Enter
最後に念のためpwdコマンドでカレントディレクトリを確認します。
foo[bar]% pwd Enter
/home/appi/foo/integral
foo[bar]%
4.3. ソースの雛型のコピー
求積問題でも雛型を用意してあります。
以前と同様に、次のようにcpコマンドを用いて、ファイルを自分のホームディレクトリにコピーします。
foo[bar]% cp ~yoshiko/skel/integral/* . Enter
foo[bar]%
ファイルをコピーしたら、念のためlsコマンドで確認してみてください。
4.4. ソースの編集
乱数の度数分布で起動したMuleが、まだ動いている場合はそのまま使用しても構いません。終了してしまった場合はもう一度Muleを起動してください。
ファイルcrude.c
とhitormiss.c
の二つを編集して、プログラムを完成させます。
また被積分関数については、ファイルFunction.c
を編集して定義します。
4.5. 実行ファイルの構築
求積問題でもMakefileを用意してあります。ファイルcrude.c
とhitormiss.c
をコンパイルするには、ktermで次のように入力します。
foo[bar]% make Enter
プログラムに間違いがなければ次のように出力され、自動的に実行ファイルcrude
とhitormiss
が構築されます。
foo[bar]% make
ln -s ../rand-distribution/Summation.h .
ln -s ../buffon/Trial.h .
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c crude.c
ln -s ../rand-distribution/Summation.c .
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c Summation.c
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c Function.c
gcc -o crude crude.o Summation.o Function.o -lm
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c hitormiss.c
ln -s ../buffon/Trial.c .
gcc -Wall -O2 -pipe -D__USE_FIXED_PROTOTYPES__ -c Trial.c
gcc -o hitormiss hitormiss.o Trial.o Function.o -lm
ln -s ../rand-distribution/rand .
ln -s ../rand-distribution/graph.sh .
foo[bar]%
エラーや警告が出た場合は、内容を調べて修正し、もう一度makeしてみます。
正常に構築できたようであれば、次に進みます。
4.6. プログラムの実行テスト
4.6.1. crudeコマンド
実行する前に、パラメータ
foo[bar]% ./rand 10 | ./crude 0 1 Enter
次のように10ステップ実行されます(値は異なるかもしれません)。
foo[bar]% ./rand 10 | ./crude 0 1
1 0.857868
2 0.921152
3 0.931160
4 0.909657
5 0.791601
6 0.823858
7 0.807871
8 0.828640
9 0.833127
10 0.849034
foo[bar]%
正常に動作しなかったり、値がめちゃくちゃだったりした場合は、Muleで修正して再びmakeします。
4.6.2. hitormissコマンド
実行する前に、パラメータ
foo[bar]% ./rand 10 | ./hitormiss 0 1 1 Enter
randコマンドは10個の乱数を発生させますが、hitormissコマンドは1ステップあたりに2つの乱数をとるので、次のように5ステップだけ実行されます(値は異なるかもしれません)。
foo[bar]% ./rand 10 | ./hitormiss 0 1 1
1 1.000000
2 1.000000
3 1.000000
4 1.000000
5 1.000000
foo[bar]%
正常に動作しなかったり、値がめちゃくちゃだったりした場合は、Muleで修正して再びmakeします。
4.7. グラフの作成
試行回数graph.sh
が用意してあります。
4.7.1. crudeコマンド
試行回数が10000回で、それぞれ
foo[bar]% ./rand 10000 | ./crude 0 1 | sh graph.sh -0 crude Enter
また、試行回数10000回目での積分値を数値で知りたい場合は、次のように入力することで知ることができます。
foo[bar]% ./rand 10000 | ./crude 0 1 | tail Enter
9991 0.784151
9992 0.784167
9993 0.784184
9994 0.784204
9995 0.784208
9996 0.784210
9997 0.784216
9998 0.784232
9999 0.784196
10000 0.784164
foo[bar]%
入力するとxgraphのウィンドウが現れますので、収束している様子を確認したら、印刷して終了します。
4.7.2. hitormissコマンド
試行回数が10000回で、それぞれ
foo[bar]% ./rand 20000 | ./hitormiss 0 1 1 | sh graph.sh -0 hit-or-miss Enter
randコマンドに与える引数は、試行回数の2倍でなければならないことに注意してください。入力するとxgraphのウィンドウが現れますので、収束している様子を確認したら、印刷して終了します。
また、試行回数10000回目での積分値を数値で知りたい場合は、次のように入力することで知ることができます。
foo[bar]% ./rand 20000 | ./hitormiss 0 1 1 | tail Enter
9991 0.783906
9992 0.783927
9993 0.783949
9994 0.783970
9995 0.783992
9996 0.784014
9997 0.784035
9998 0.784057
9999 0.783978
10000 0.783900
foo[bar]%