プログラム合宿2009春報告

概要

2009年3月16日(月)から19日(木)にかけて行ったプログラム合宿2009春の内容のまとめ兼報告です。

1日目 03/16(月)

↑top
↓

1日目はまず、それぞれのPCでCプログラミングができるよう環境構築し、コンパイルの仕方や基本文法などの基礎的なことを丁寧に扱いながら、実際に微分積分の数値計算をしてもらいました。

環境構築

ほぼ全員の参加者がCコンパイラを持っていなかったため、こちらが予め用意しておいたインストーラを配布して環境構築作業を行いました。コンパイラはGCCで、Equation Solutionに置いてあるwin用インストーラの20090206版(gcc-4.4-20090206-32.exe)を利用しました。

次に、プログラム作成から実行まで。流石にテキストエディタの説明は不要でしたが、コンパイル、プログラムの実行の仕方については、やはり存在感がなく使いづらいWindowsのコマンドプロンプトについていくらか説明が必要でした。

Windows標準のコマンドプロンプトはcmd.exeという名前で、"cmd"とファイル名を指定して実行することで起動できます。Vistaなら、フォルダのアドレスバーで"cmd"と入力すれば、そのフォルダが作業フォルダになった状態でコマンドプロンプトが立ち上がるので便利です。

基礎文法

C言語の基礎文法について簡単に扱いました。ここでいちいちその内容について詳しくは書きません。

テンプレート(おまじない)
1. #include
2. int main(void){～}
関数(型B function(型A a, ...){～})
変数(型A x=初期値;)
出力(printf命令)

実際には、次の微分積分の数値計算のプログラムを実際に組む演習と並行して必要な時点で解説しました。

数値計算(微分積分)

数値計算とは、数学的・物理的な厳密解を求める(解析的に解く)ことではなく、近い値の解を求める(数値的に解く)ことを言います。特に現代では、近似解をコンピュータを用いて計算することを指すことが多いです。

微分積分の数値計算においては、極限の考え方で計算します。具体的には、次のような式で近似します。

微分係数: $f'(a)≒{f(a+h)-f(a)}/h (h<<1)$
定積分: $∫[x=a~b]f(x)dx≒∑[k=1~n]f(x_k)Δx (Δx=(b-a)/n<<1, x_k=a+kΔx)$

このように、hやΔxが十分小となるように取ることで近似する考え方で、三角関数の微分や積分をしてもらいました。

困ったことに、刻み幅が大きすぎると近似精度が悪くなり、逆に小さすぎてもコンピュータの計算許容範囲を越えてしまいます。こういった問題点を解消するため、より精度のよいアルゴリズムや近似式がいくつかありますが、単純のためこれ以上は扱いませんでした。

2日目は1日目の発展として、さらに実用的な微分方程式の数値的解法を扱いました。

オイラー法

オイラー法とは、(1階の常)微分方程式とその初期条件が与えられた時に、発展的にその数値解を求める最もシンプルな方法のひとつです。(単純ゆえに精度は非常に悪い)

例えばある物体(位置x)の速度が分かっている、つまり微分方程式 $dx/dt=f(x,t)$ が与えられている時を考えます。初期条件(時刻t=t_{_0}における位置x=x_{_0})が与えられていれば、その物体の位置の時間発展を次のようにシミュレートすることができます。微分の近似式を等式変形しただけですが、時刻tにおける物体の位置xと速度f(x,t)から、時間Δt後の物体のおよその位置が分かります。

$x(t+Δt)=x(t)+Δt*f(x,t) (Δt<<1)$

1日目同様に三角関数を例として計算してもらいました。当然、同じ周期を保って振動する結果が得られるのですが、1次近似のために精度が悪いということも確認してもらいました。

見える化(Visualize)

せっかくなので、結果をCSV出力してエクセルのグラフにする方法も扱いました。

while ( t < 10.0 ) {
  printf( "%f,%f", t, x );
  x = x + dt * f(x,t);
  t = t + dt;
}

例えば上のように、ループ内のprintf文でカンマ区切りに値を表示するようにしておきます。このプログラムをコンパイルして、結果を出力するには次のようにします。

*コマンドプロンプト*での入出力例
コマンド1	入力	`a.exe > output.csv`
コマンド1	出力	※コマンドプロンプト上には何も表示されず、output.csvに結果が保存される

標準ではCSVファイルはエクセルで開くよう設定されているので、出力されたCSVファイルをダブルクリックすればエクセルで開きグラフを作ることができます。

4次のルンゲ=クッタ法

今までの考え方は精度が非常に悪いものであるため、実際に利用されることはあまりありません。それに対して、(4次の)ルンゲ=クッタ法(RK4)と呼ばれる方法は精度がよく、この方法で(常)微分方程式を解くと解の発散をより抑えることができます。

2日目の時間が残り少なかったこともあり、当日の合宿参加者でRK4を実装させることができたのは1人だけでした。

3日目 03/18(水)

↑top
↑
↓

3日目は、乱数を用いて大数の法則とモンテカルロ法を扱いました。

乱数

C言語にはrand()関数が標準で用意されていて、0から32767(16進数で0x7FFF)までの範囲の整数をランダムに返します。[0,1)区間の乱数がこれから必要なので、まずは[0,1)区間の実数をランダムに返す関数を作り、使用してもらいました。

double frand(void) {
  return (double)rand()/0x8000;
}

モンテカルロ積分

大数の法則により、試行を何回も繰り返して事象が起きた回数を総試行回数で割れば、事象の確率に近づきます。これを利用して積分の値を求めるのがモンテカルロ積分です。合宿では、モンテカルロ法により単位四分円の面積からπの大雑把な値を求めました。

[0,1)^{^2}の元Pをランダム取る試行を考えます。このとき、点Pが単位四分円D := {(x,y) | x,y>0, x^{^2}+y^{^2}≦1} に含まれる事象の確率は、面積比と等しいπ/4になります。よって、試行を十分な数だけ行って確率を計測し、4倍すればπの大雑把な値を求めることができます。

このように、モンテカルロ法を用いて(重)積分の値を比較的簡単に求めることができます。ただし、一様な乱数で十分な試行回数を行わなければ正確な値は得られず、1次近似ほど効率はよくありません。

比較的全員すんなり完成させられたので、試行回数をさらに増やして計算してもらったり、四分円内部の点Pの座標をCSV出力してエクセルで散布図を描くなどしてもらいました。

気体分子運動モデル

モンテカルロ法を用いれば、確率モデルの時間発展を調べることができます。例えば次のような気体分子運動モデルを考えます。

二つの容器V_{_1},V_{_2}(同体積)があり、単位時間当たり1個の気体分子しか通過できない通路(穴)で繋がれているとする
気体分子は全てでN個あり、V_{_n}の分子数をN_{_n}とする
初めはN_{_1}=N, N_{_2}=0
単位時間当たりに気体分子がV_{_1}からV_{_2}へ移動する確率 $p_1=N_1/N$

frand()で得た乱数がp_{_1}以下であればV_{_1}からV_{_2}へ、そうでなければV_{_2}からV_{_1}へ1個の分子が移動したと考えます。このような試行を十分数繰り返すことで、気体分子の運動の時間発展、つまり気圧平衡をシミュレートすることができます。

例としては単純化したモデルを扱いましたが、より厳密なモデルを採用してももちろん同じように調べることができます。あまり時間はありませんでしたが、容器を3つにするなどの他のモデルでも試してもらいました。

4日目 03/19(木)

↑top
↑
↓

4日目はセルオートマトン・CGを扱いました。本当はCやRubyでCGを扱えるとよかったのですが、ライブラリを用意したりする時間がなかったので、Schemeメインです。

セルオートマトン

簡単に言うと、セルオートマトンはライフゲームのより一般的・学術的な総称です。ライフゲームは普通2次元ですが、今回は1次元のセルオートマトンを扱いました。

一列にセルが並んでいて、個々のセルは生存・死亡のどちらかの状態を取ります。ある時刻において「生存」のセルは、もし両隣のセルが「生存」ならば次の時刻には「死亡」する(過密死)というようなルールを作ります。そしてセルの初期状態を与え、ルールに基づいた時間発展を図示します。するとルールによってはフラクタル図形が出来上がったり、色々と興味深い結果を得られます。

このルールは、セルの生存を1、死亡を0と同一視(二進表示)することで、簡単に配列で定義することができます。ルールの例を下の表に示します。この表では、生存セルを「■」、死亡セルを「□」で図示しています。

発展ルールと二進表示
インデクス	二進表示	セル	次セル
0	`0b000`	□□□	□
1	`0b001`	□□■	■
2	`0b010`	□■□	■
3	`0b011`	□■■	□
4	`0b100`	■□□	■
5	`0b101`	■□■	■
6	`0b110`	■■□	□
7	`0b111`	■■■	□

図1: セルオートマトンの時間発展図

このルール例の場合、{0,1,1,0,1,1,0,0}という配列(Schemeの場合なら'(0 1 1 0 1 1 0 0)というリスト)を作ればいいのです。現在のセル及び隣接するセル(3列目)の情報を二進表示(2列目)と見なし10進数の整数に直し(1列目)、ルール配列(リスト)のその番号の要素を次のセルの状態(4列目)として更新すればいいのです。

もちろん、C言語などで標準出力することで時間発展を調べることもできます。しかしどうせならグラフィックスにしてみたいですよね。（図1は、上のルールで1個のセルだけ生きている初期状態から始めたもの）

CG (with Scheme)

基礎情報処理演習でSchemeに慣れている参加者が多かったので、基礎は省いてSchemeでグラフィックスを扱う方法の復習から入りました。Schemeでグラフィックスを扱うためには、まずプログラムの始めに次のように書きます。

(require (lib "graphics.ss" "graphics"))
(open-graphics)
(define w (open-viewport "Cellular Automaton" 500 500))

そして、wに対して線, 四角, 楕円, 多角形を描くには、次のような命令を使います。なお図形命令で内部を塗り潰す場合は、draw-～ではなくdraw-solid-～という命令を使います。

((draw-line w) (make-posn x1 y1) (make-posn x2 y2) "color")
((draw-rectangle w) (make-posn x y) width height "color")
((draw-ellipse w) (make-posn x y) width height "color")
((draw-polygon w) (list (make-posn x1 y1)
                        (make-posn x2 y2)
                         ;
                        (make-posn xn yn) )
                  (make-posn x0 y0) 
                  "color" )

合宿では線の描き方しか教えませんでしたが……f(^^;

番外編１～課外授業～

↑top
↑
↓

合宿1日目が終了した後、2日目に参加できない人のために課外授業を行いました。2人の参加者がRK4を実装させるところまで行きました。

また、Schemeで数値計算をする際気をつけると得する仕様も教えました。Schemeでは実数型と有理数型が明確に分かれていて、C同様に有効桁数が存在する実数型と違い、有理数型ではデータが2つの整数(分子/分母)に分けて扱われます。さらにSchemeでは内部的に仮想メモリを利用しており、C言語と違い莫大な桁の整数データ（但しその分遅い）を扱うことができます。ゆえに有理数型の精度も非常に巨大です。

計算中に一度でも実数型の値を用いるとその結果も実数型になりますが、有理数型のみで計算をしている限りは計算結果は有理数型であり、従って計算の精度を有理数型の精度に保つことが可能です。刻み幅を十分小さく取らなければ精度を得られない微分積分の計算において、この性質は非常に強力です。実数型で刻み幅を小さくしすぎるとすぐに値が計算不可能な範囲に達するのに対し、有理数型はその計算可能範囲が尋常ではないのです。

精度の悪い実数型の値を精度の良い有理数型へ変換するためには、inexact->exact命令を使用します。この命令は実数以外でも、不正確な値を正確な値へ変換するために使うことができます。

(inexact->exact pi)
(inexact->exact 1.25+i5.00)

番外編２～個別指導～

↑top
↑
↓

他にも合宿中、質問を受けたりしたときに以下のようなことを教えました。

[C] FILE型の扱い、標準入出力との違い
[C] 配列とポインタの扱い・考え方
[C] 型の強制変換(キャスト)
[Scheme] lambdaと関数適用
[Ruby] 添削、キーボード入力、フラグの使い方
テントマップによる乱数生成アルゴリズム

また合宿中、余裕のある参加者は扱った内容以外にも様々なプログラムを自由に(自発的に)作っていました。

[VB] ライフゲーム
[VB] シューティング
[Ruby] シューティング
[Ruby] (ごく簡単な)RPG
[Scheme] sinの級数展開のグラフ表示

資料

↑top
↑
↓

暇があったら今回の合宿で参考にした資料、作られたプログラム等もまとめる予定。

図2: 組込のlog関数と級数展開で近似したlog関数のグラフ
図3: 組込のsin関数と級数展開で近似したsin関数のグラフ

dsinx.c: [UDA] sinxのx=πにおける微分係数を求める
sum.c: [UDA] πに収束する無限和の計算
euler.c: [UDA] オイラー法でx(0)=0,v(t)=cos(t)の時間発展を調べる
kepler.c: [UDA] RK4でケプラー問題を数値的に解く
calc_pi.ss: [UDA] 有理数の四則演算のみでπを計算する
ame.c: [あずまん] モンテカルロ法でπを求める
para.ss: [@b] Schemeで媒介変数表示のグラフを描画する
dif.ss: [@b] n回微分を計算する; 図2, 図3は、dif.ssで級数展開の係数を求め、para.ssでグラフを描画している様子

まとめ

↑top
↑

1-3日目の内容は主に後期理学部科目「物理学情報処理論1」の内容を参考にしています。Windows環境を基本にした合宿と、学校のLinux環境を基本にした物理情報では内容のずれが割とありますが、今回の合宿で参加者は物理情報の単位取得に必要な最低限のスキルは身につけたと言えるでしょう。また4日目の内容も、敢えて授業と関連付けるならば前期一般教養科目「コンピュータグラフィックス実習A」で同様の内容を扱っています。まぁCG実習はSchemeでなくRubyを使用しているので、かなり勝手が異なりますが。

今回は基本中の基本を中心に扱いましたが、参加者にはかなりプログラミング力をつけてもらえたのではないかと思っています。新年度のS2Sの方向性にもよるのでまだ確定はしていませんが、夏もまた合宿をやると思われます（こちらは本当に合宿になるかも！？）。今後ともS2Sをよろしくお願いします。

S2S