プログラミング演習Ⅰ 第 14 回

2015/07/10
プログラミング演習Ⅰ
第 14 回
総合演習(2)
今日の目標:今までの命令文を使いこなしてまとまったプログラムを書く。
1. 乱数を使った数値積分
今週は乱数を用いて関数積分を行う。これはモンテカルロシミュレーション法と呼ばれる。2次元
の関数の積分の場合、例えば[0,1]に分布する乱数を 2 組用意し、各々を(x, y)の組として扱う。
理想的には図1に示す象現内に均等に分布するはずである。そこで、この象現内である関数 f(x)
より大きい組の数と全数の比は、f(x)の積分値(f(x)より下にある面積)と象現全体の面積(こ
の場合は x, y 共に[0,1]なので、1に等しい)の比になるはずである。よって f(x)の積分値

1
0
f ( x) dx
1
0.9
が求まる。
0.8
0.7
0.6
0.5
まず図1に示す
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
プログラム例:
図 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
2015/07/10
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) プログラム例1を入力、コンパイルし、動作を確認せよ。
(2) (1)のプログラムを改良して、図
2に示すように円の 1/4 の面積
を求め、これよりを計算するプ
ログラムを作成せよ。この時、
1000 組数毎に中間結果を出力し、
収束の度合いを見よ。
(ヒント)1/4円内では
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
cos x dx
1000 組数毎に中間結果を出力し、収束の度合いを見よ。この積分の理論値と比較せよ。
注意 1:cos を使用するので,#include <math.h>を入れておくことと、コンパイル時
プログラミング演習第 14 回
2
2015/07/10
-lm を含めることを忘れずに。
注意 2:積分範囲が 0  x 
の面積は

となる。
2

2
であることに注意(y は 0 
y  1 で不変)。よって象現全体
(5) (1)のプログラムを改良して、以下の積分値を求めよ。
1 x
0 e
dx
1000 組数毎に中間結果を出力し、収束の度合いを見よ。この積分の理論値と比較せよ。
注意:e x を使用するので,#include <math.h>を入れておくことと、コンパイル時-lm
を含めることを忘れずに。
z
(6) (1)のプログラムを改良して、半径 1.0
の球の体積を求めなさい。
1.0
1000 組数毎に中間結果を出力し、収束の
度合いを見よ。この積分の理論値と比較
せよ。(ヒント)体積を求める場合は 3
組
の
乱
数
x,
y,
z
( 0  x  1, 0  y  1, 0  z  1)を発生
x
1.0
する。一辺が 1.0 の直方体の中の 1/8
球 を 考 え る 。 こ の 球 内 で は
x 2  y 2  z 2  1 が成り立つ。直方体内の
1.0
全点数の内球内の点数を求めて、体積を
求めよ。
y
(7) (1)のプログラムを改良して、高さと底
z
面の半径が 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 軸上にあるとする。こ
1.0
の円錐内では x  y  (1  z ) が成り立
y
2
2
x
2
つ。直方体内の全点数の内円錐内の点数
を求めて、体積を求めよ。
(8) (発展)余裕のある人はやってみてくだ
さい。
πはライプニッツの公式でも計算出来る。

1 1 1 1
 1   
4
3 5 7 9
(1) n

n  0 2n  1

(2)でモンテカルロ法で求めたπとライプニッツ公式で求めたπ、およびその絶対誤差
(|推定値-真値|)を表示し、比較せよ。πの真値として M_PI を用いよ。(2)と同様、1000
組数毎に中間結果を出力し、収束の度合いを見よ。
プログラミング演習第 14 回
3
2015/07/10
出力例:
$ ./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/15)16 時までとします。電情系事務室(7 号館 2 階 221 号室)内の
プログラミング演習用メールボックスに提出せよ。
注意:提出レポートで明らかに動作しないプログラムは再提出してもらうので、十分デバッグ、
動作確認すること。再提出、期限外はその度に評価を1ランク下げるので、注意。また、未提出、
ならびに再提出を未提出のものは 0 点とする。学生番号、氏名、出題日、課題番号を各プログラム
第 1 ページ先頭にコメント行を使って明記のこと。複数枚にわたる場合はホチキス等でとめること。
簡単な課題なので、必ず自分で演習に取り組むこと。
プログラムにはコメントを十分記入し、インデント、改行を適当に入れて読みやすいプログラム記入
を心掛けること。
動作しないプログラム、指示どおりの機能しないプログラムの他、適切なコメントやインデント
のない提出プログラムも再提出の対象となるので注意。
前回再提出となった課題の期限も来週水曜(7/15)の 16 時とします。再提出とされた元のレポー
トの上に修正したレポートを一緒に綴じて提出のこと。今週のレポートと再提出のレポートは一緒
に閉じないこと。
今回のレポートは7月17日、5・6校時に 7-228 において返却予定。このレポートの再提出
は 7 月 22 日の 16 時電情系事務室(7 号館 2 階 221 号室)内のプログラミング演習用メールボッ
クスとする。これ以降は一切のレポートは受け付けない。
プログラミング演習第 14 回
4