2014/07/18
プログラミング演習Ⅰ
第 14 回
総合演習（2）
今日の目標：今までの命令文を使いこなしてまとまったプログラムを書く。
1. 乱数を使った数値積分
今週は乱数を用いて関数積分を行う。これはモンテカルロシミュレーション法と呼ばれる。２次元
の関数の積分の場合､例えば[0,1]に分布する乱数を 2 組用意し､各々を(x, y)の組として扱う。
理想的には図１に示す象現内に均等に分布するはずである。そこで、この象現内である関数 f(x)
より大きい組の数と全数の比は、f(x)の積分値（f(x)より下にある面積）と象現全体の面積（こ
の場合は x, y 共に[0,1]なので、１に等しい）の比になるはずである。よって f(x)の積分値

1
0
f ( x) dx
1
0.9
が求まる。
0.8
0.7
0.6
0.5
まず図１に示す
0.4
1
 x dx
0
0.3
( f (x) = x)
0.2
0.1
を求めてみる。
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
プログラム例：
図１
/**************************************************
モンテカルロ法による x=y の積分
***************************************************/
#include <stdlib.h>
#include <stdio.h>
#define NUMPOINTS 10000
float unitrand(void); /* [0, 1]の乱数発生 */
int main(void)
{
int i, count;
float x, y;
count = 0;
for (i = 0; i < NUMPOINTS; i++) {
x = unitrand();/* [0, 1]の乱数一組(x, y)発生 */
y = unitrand();
/* 発生乱数対が積分領域内ならカウンタインクリメント */
if (x > y) /* f(x) > y ? */
count++;
}
プログラミング演習第 14 回
1
2014/07/18
printf("x = y の面積は
%f¥n", (float)count/NUMPOINTS);
return 0;
}
/* [0,1]の乱数発生 */
float unitrand()
{
return (float)rand()/RAND_MAX;
}
出力例：
$ ./example14_1.exe[Enter]
x = y の面積は 0.499300
2. 本日の演習
(1) プログラム例１を入力、コンパイルし、動作を確認せよ。
(2) (1)のプログラムを改良して、図
２に示すように円の 1/4 の面積
を求め､これよりを計算するプ
ログラムを作成せよ。この時、
1000 組数毎に中間結果を出力し､
収束の度合いを見よ。
（ヒント）１／４円内では
x 2  y 2  1 が成り立つ。
1
0.9
0.8
0.7
0.6
0.5
出力例：
$ ./lab14_2.exe[Enter]
0.4
1000, πの推定値は 3.160000
0.3
2000, πの推定値は 3.136000
3000, πの推定値は 3.149333
0.2
4000, πの推定値は 3.152000
0.1
5000, πの推定値は 3.167200
6000, πの推定値は 3.170000
0
0
0.1
0.2
0.3
7000, πの推定値は 3.168571
8000, πの推定値は 3.169500
9000, πの推定値は 3.176444
πの最終推定値は 3.178800
(3) (1)のプログラムを改良して、以下の積分値を求めよ。
1 2
x
0

0.4
0.5
0.6
0.7
0.8
0.9
1
図 2
dx
1000 組数毎に中間結果を出力し､収束の度合いを見よ。この積分の理論値と比較せよ。
(4) (1)のプログラムを改良して、以下の積分値を求めよ。


2
0
sin x dx
1000 組数毎に中間結果を出力し､収束の度合いを見よ。この積分の理論値と比較せよ。
注意 1:sin を使用するので,#include <math.h>を入れておくことと、コンパイル時
プログラミング演習第 14 回
2
2014/07/18
-lm を含めることを忘れずに。
注意 2:積分範囲が 0  x 
の面積は

となる。
2

2
であることに注意（y は 0 
y  1 で不変）。よって象現全体
(5) (1)のプログラムを改良して、以下の積分値を求めよ。
1 x
0 e
dx
1000 組数毎に中間結果を出力し､収束の度合いを見よ。この積分の理論値と比較せよ。
注意：e x を使用するので,#include <math.h>を入れておくことと、コンパイル時-lm
を含めることを忘れずに。
(6) (1)のプログラムを改良して、半径 1.0
ｚ
の球の体積を求めなさい。
1000 組数毎に中間結果を出力し､収束の
1.0
度合いを見よ。この積分の理論値と比較
せよ。（ヒント）体積を求める場合は 3
組
の
乱
数
x,
y,
z
（ 0  x  1, 0  y  1, 0  z  1）を発生
する。一辺が 1.0 の直方体の中の 1/8
球を考える。この球内では
1.0
x
x 2  y 2  z 2  1 が成り立つ。直方体内の
全点数の内球内の点数を求めて、体積を
求めよ。
1.0
y
(7) (1)のプログラムを改良して、高さと底面
ｚ
の半径が 1.0 の円錐の体積を求めなさい。
1000 組数毎に中間結果を出力し､収束の
度合いを見よ。この積分の理論値と比較
せよ。
（ヒント）体積を求める場合は 3 組
の
乱
数
x,
y,
z
（ 0  x  1, 0  y  1, 0  z  1）を発生
1.0
1.0
する。一辺が 1.0 の直方体の中に円錐の
1/4 が入っているとする。底面は xy 平面
にあり、頂点は z 軸上にあるとする。こ
の円錐内では x 2  y 2  (1  z ) 2 が成り立
つ。直方体内の全点数の内円錐内の点数
を求めて、体積を求めよ。
(8) （発展）余裕のある人はやってみてくだ
さい。
πはライプニッツの公式でも計算出来る。

1 1 1 1
 1   
4
3 5 7 9
x
1.0
y
(1) n
n  0 2n  1


(2)でモンテカルロ法で求めたπとライプニッツ公式で求めたπ、およびその絶対誤差
(|推定値－真値|)を表示し、比較せよ。πの真値として M_PI を用いよ。(2)と同様、1000
プログラミング演習第 14 回
3
2014/07/18
組数毎に中間結果を出力し､収束の度合いを見よ。
出力例：
$ ./lab14_8.exe[Enter]
πの推定値(1000 点)：3.136000 誤差 0.005593 (MC) 3.142592 誤差 0.000999
(Leibniz)
πの推定値(2000 点)：3.154000 誤差 0.012407 (MC) 3.142090 誤差 0.000497
(Leibniz)
（中略）
πの最終推定値：
3.xxxxxx 誤差 0.xxxxxx (MC) 3.xxxxxx 誤差 0.xxxxxx
(Leibniz)
(2)、（4）
、(6)の各リストと、実行時の表示印刷(script 使用)を提出せよ。理論値（自分で
積分してみよ）も計算とともに示しなさい。できた人は(8)も表示印刷ともに提出せよ。
提出期限は来週水曜（7/23）16 時までとします。電情系事務室（7 号館 2 階 221 号室）外の
プログラミング演習用メールボックスに提出せよ。
注意：提出レポートで明らかに動作しないプログラムは再提出してもらうので、十分デバッグ、
動作確認すること。再提出､期限外はその度に評価を１ランク下げるので、注意。また、未提出、
ならびに再提出を未提出のものは 0 点とする。学生番号、氏名、出題日、課題番号を各プログラム
第 1 ページ先頭にコメント行を使って明記のこと｡複数枚にわたる場合はホチキス等でとめること。
簡単な課題なので､必ず自分で演習に取り組むこと。
プログラムにはコメントを十分記入し､インデント、改行を適当に入れて読みやすいプログラム記入
を心掛けること。
動作しないプログラム､指示どおりの機能しないプログラムの他、適切なコメントやインデント
のない提出プログラムも再提出の対象となるので注意。
前回再提出となった課題の期限も来週水曜（7/23）の 16 時とします。再提出とされた元のレポー
トの上に修正したレポートを一緒に綴じて提出のこと。今週のレポートと再提出のレポートは一緒
に閉じないこと。
今回のレポートは期末試験において返却予定。このレポートの再提出は 7 月 30 日の 16 時電情
系事務室（7 号館 2 階 221 号室）外のプログラミング演習用メールボックスとする。これ以降は一
切のレポートは受け付けない。
プログラミング演習第 14 回
4