情報 課題
4
田浦健次朗
2010.6
月某日
1
急がば回れ
:
間違い探し
課題を始める前にいくつか「よくある間違い」を事前に経験してもらうためのプログラムを, 授業 HP に掲載している. 急がば回れで実際の課題に取り組む前に以下を行うこと. • 講義 HP で「課題 4 で用いる題材 (間違ったプログラム集, PGM/PPM 画像サンプル)」のリン クをつついてダウンロードページを見る. • そこから各プログラムをダウンロード • エディタで開き, コメントを読む • python で実行し, エラーが出ることを確認 • 指示された修正を施すか, または自分で修正方法を考え, 修正・実行する2
課題実行・提出方法
以下の課題は, • kadai4 というフォルダを作り, その中にファイルを作っていく. • 基本課題, オプション課題, 考察課題とあるうち, – 基本課題については必須 – オプション課題については余力があれば実行 – 考察課題については余力があれば考察 授業時間内にも関わらず「余力がありません (のでオプション課題はやりません)」というの は無し. • 基本課題については, 節 (問題) 番号 N に対応して, exN.py というファイル名でプログラムを 作る (例えば 3 節の問題に対しては ex3.py というファイルを作る). N は 3 から始まる. • オプション課題も同様. やらなかった場合は, (当然) そのファイルを含めない. • 考察課題については別途, 節 (問題) 番号 N に対応して, 考察を記述したテキストファイル (exN.txt) を作る. 中身の形式は自由.• その他課題ごとに指示があれば指示にしたがってファイルを含める (出力結果など) • 最後に kadai4 フォルダをファインダで右クリック → 圧縮して, kadai4.zip ファイルを作って 提出する. • 授業時間内にベストを尽くせば大体の人が終わるよう, 時間を設定する予定.
3
簡単な繰り返し
基本課題
cos x の値を [−π, π] の範囲で表示する Python プログラムを書け. この区間内に適当な数 (例 えば 100 個) の点をとり, 各点における x と cos x の値を表示する (一行につき二つの数字). 以下のような概形の (一旦関数を書いて, それを呼び出す) プログラムとして書く事. def print cos():. . . . . . print cos() 実行すると以下のような出力が得られれば正解. -3.14159265359 -1.0 -3.13530946828 -0.999980260856 -3.12902628298 -0.999921044204 -3.12274309767 -0.999822352381 . . . ヒント: • for 文を使う (課題 4 のプリント 10 節) • 一行に二つの値を print するやり方 (課題 4 のプリント 8 節) • cos 関数を使うときのおまじない (課題 4 のプリント 3 節) などに気をつけながら書いてみよ. 上記の結果が確認できたら, 出力をコマンドラインのリダイレクト機能 (>) を使ってファイルに 保存する.
python ex1.py > cos.dat
このように各行に, 二つの値が並んだデータは, gnuplot というコマンドで可視化できる. gnuplot
と端末に入力して gnuplot を起動すると, 以下のプロンプト (gnuplot>) が現れる. G N U P L O T
Version 4.2 patchlevel 4 . . .
Send bug reports and suggestions to <http://sourceforge.net/projects/gnuplot>
Terminal type set to ’wxt’ gnuplot>
このプロンプト (gnuplot> ) に向かって, plot ”cos.dat”
と入力するとグラフが現れる. 上記の代わりに plot ”cos.dat” with lines
と入力すると, 線で結んだ表示になる. gnuplot を終了するには gnuplot> プロンプトに向かって, Ctrl-d (Ctrl キーを押しながら d) を入力する. 何度もデータを変えてはグラフを表示させる場合, 毎 回 gnuplot を起動・終了させる必要はない. gnuplot を開いておいて, 必要な plot コマンドを打てば グラフを再描画してくれる. gnuplot 用に一つコマンドライン端末の窓を開いておき, そこで gnuplot を立ち上げ放しにしておけばよい. または, 以下のような内容のテキストファイル (例: show.gpl) plot ”cos.dat” pause -1 を用意しておけば, 端末から以下のコマンド gnuplot show.gpl を実行するだけですぐにグラフが表示される. 端末でリターンキーを押すと終了する. グラフの表示は Excel を使ってもできた. その場合, 二つの数字の間にカンマ (,) を出力し, 出力 を.csv で終わる名前のファイルに保存すればよい (print x, ”,”, cos(x)). ただし, いろいろなデータ を生成してはそれをグラフ表示する, ということを繰り返す場合, 慣れてしまえば gnuplot (特に後 者の方法) の方が早い. その上, Excel には表示できるデータの数に割と小さな制限がある.
4
フィボナッチ数列
基本課題
n を入力とし, 次の規則で決まる数列 (フィボナッチ数列): u0 = 1 u1 = 1 un = un−1+ un−2 (n≥ 2) の, unを求める関数 fib を書け. それを呼び出して u10を表示せよ. プログラム全体の雰囲気は以下のようになるだろう. def fib(n): . . . return . . . print fib(10) ヒント: • for 文 • 変数と代入文 (課題 4 プリント 7 節) を使う いくつ変数を用意して, どう更新するかがポイント よくわからなかったら, 以下のような表をどう埋めていくか考えるとよい. 手で埋めるだけでも 何をしたらよいか分かるかもしれないし, わからなければ Excel の数式機能を使って実際に求めて みてもよい. それができればプログラムに直すのはおそらく簡単である. i ui−2 ui−1 ui 2 1 1 3 .. . n5
基本課題
上のプログラムを少し手直しして, n が与えられたら, u2, u3, . . . , un までを順に表示するプロ グラムに変更せよ. 各行の第一列に i, 右側に uiを表示せよ.ヒント:
• 上記のプログラムを少し手直しして, for 文を実行しながら途中経過を print する文をはさめば
良い
プログラムの概形はおそらく以下のようになり, print fib の中で print 文が呼ばれる. def print fib(n):
. . . print fib(10) そのプログラムを実行すると以下のような出力が得られるはずである. 2 2 3 3 4 5 5 8 6 13 7 21 8 34 9 55 10 89
6
一般の漸化式
せっかく作ったプログラムにフィボナッチ数列だけを計算させるのはもったいない. 少し書き換え ればより一般的な漸化式の計算ができるようになる. 関数を「パラメータ化」して用途を広げるの は常套手段.オプション課題
係数 a, b, 最初の 2 項 u0, u1と n ≥ 2 を受け取り, 以下の漸化式で定まる数列の n 項目 (un) ま でを, 前問と同様の形式で表示する関数を書け. もちろん, a = b = −1, u0 = u1 = 1 とすればそれ はフィボナッチ数列となる. un+ aun−1+ bun−2 = 0 (n≥ 2) 予想される概形は以下.def print recurrence(a, b, u0, u1, n): . . .
. . .
7
前問の数列に関して, • a = 1, b = 1 とすると, この数列は 3 を周期とした繰り返しになる. • a = 0, b = 1 とすると, 周期は 4 になり, • a = −1, b = 1 とすると, 周期は 6 になる. これは実際にやって見れば観測できる.考察課題
さて, 周期を 5 とするには a, b をどのように選んだらよいか? それを求め, 実際に確かめてみ よ (gnuplot で可視化してみると感じが分かる). 周期 13 などという器用なこともできるのか? ヒント この数列の一般項は, 漸化式の特性方程式 x2+ ax + b = 0 の 2 解 α, β によって決まる. α = β の場合は若干面倒なので省略するが, α6= β ならば, un= Aαn+ Bβn (A, B は u0, u1に応じて定まる定数) となる.8
微分方程式
(漸化式を解くのと, 微分方程式を解くのは似ている).オプション課題
係数 a, b, 初期値 u(0), u0(0), および t を受け取り, u(t) の近似値を計算する関数を書け. u00(t) + au0(t) + bu(t) = 0 予想される概形.def solve de(a, b, u0, u 0, t): . . . solve de(0, 1, 1.0, 1.0, 10.0) 自分で適当な刻み幅 ∆t を設定し, 0 から t における各地点での値を, 前問と同様の形式 (第一列 に t, 第二列に u(t)) で表示し, それを gnuplot で表示してみよ. 最後の行で a, b に渡す値を変えて (つ まり, a, b を変えて), 解の形がどのように変わるかを観測し, 自分の理解 (または直感) とあっている か確認してみよ.
ヒント: • Excel でやったときの計算の原理を思い出す. – ある t に対して u(t), u0(t) がわかれば, u(t + ∆t) が近似計算できる – ある t に対して u0(t), u00(t) がわかれば, u0(t + ∆t) が近似計算できる – 一方ある t に対して u0(t), u0(t) がわかれば微分方程式により u00(t) がわかる – 以上を合わせると, ある t に対して u(t), u0(t) がわかれば, u(t + ∆t), u0(t + ∆t) が近似計 算できる. これを繰り返して u(0), u0(0) から先の u(·) をどんどん近似計算していける. • Excel では上記を表の形式に置き, 時刻が進んだときの結果を表の下の方に追記していったが, プログラミング言語では変数を用いて上記の値を計算していく. 背景 (1) ここで述べた微分方程式の数値解法は非常に一般的なもので, 式の形から階が閉じた式と して (手動で) 計算できない場合でも通用する. つまり, u00(t) + u0(t)2+ log u(t) = 0 であろうが,
u000(t) + sin u(t) + log t = 0
であろうが, ほぼ同じやり方で解の近似値が求まってしまう. 複数の式が絡み合った方程式系でも 解法はほとんど同じであるし, u が t のみの関数ではなく, u(x, y, t) のように複数の変数の関数の時 (偏微分方程式) も, 似たやり方が適用できる. 物理, 化学, 機械工学, 電気工学, . . . 多くの分野で多か れ少なかれ, 目にすることになるだろう. 背景 (2) 特にこの形の微分方程式は t を時刻, u(t) を時刻 t における物体の位置としたときの運動 方程式や, LCR を含んだ回路の方程式として良く現れる. • a = 0, b > 0 なら単振動 (ma = −kx) ⇒ 同じ振幅を保ったまま振動する. • a > 0 ならば, a は速度に比例した抵抗を決める (ma = −av − kx) ⇒ 減衰しながら振動する (振動しない場合もある). b > 0 ならば, 直感的には b はバネ定数ということであり, 絶対値が大きければ, 多少の減衰はものと もせずにしばらく振動を続けるはずである. さらさらした水の中でバネにつながれた重りがゆらゆ らしている状態を想像すればよい. a > 0 でその絶対値が大きければ, 抵抗が大きいということで, 振動することもなく止まってしまう. どろどろ・ネバネバした液体の中でバネにつながれた重りが ゆらゆらしている状態. 一般にはこの微分方程式の解は, 特性方程式: x2+ ax + b = 0 の 2 解 α, β の値に応じて決まる. ここでも α = β となる場合は話がやや面倒になるので省略するが, α6= β であれば, この解は u(t) = Aeαt+ Beβt となる (A, B は u(0), u0(0) に応じて決まる定数). α, β が実数でない場合もこの公式は正しい (授業 で述べたとおりその場合, eαtは指数関数と三角関数の積になる). 発散する, 単振動する, 減衰振動する, 振動しないで減衰, などの挙動は α, β の値でクリアに分類 できる.
-250 -200 -150 -100 -50 0 50 100 150 200 250 0 100 200 300 400 500 600 700 800 900 1000 "b.dat"
9
前問の微分方程式に関して, 余力があったら以下を考えてみよう. • 微分方程式のシミュレーションなんて難しそう. . . • しかし, 実際やってみると意外と簡単だった という人に, 「いや, 実は意外と難しいところもある」ということも知ってもらうための話. a = 0, b > 0 という設定では運動は単振動で, 同じ振幅の運動を永遠に続けるはずである. しか し実際に長い間この計算を続けると結果は図のように, あるところまで行くと発散してしまう. 少 しやってみると分かるが, b をいろいろにかえても必ず発散する. 刻み幅 ∆t を小さくすると, 発散が 起きるまでの時間は長くなるが, それでもやがて発散する. もちろん微分方程式を近似的に解いているので, もともと厳密な解とピタリ一致する解を出せる とは期待できない. また, 一般論として近似の上に近似を積み重ねていけばどんどん誤差が増大し て, やがて現実の解とかけ離れた値になっていく, というのもある程度避けられない現象である. し かし, 振幅がフラフラしながらも一応振動は続けるとか, 逆に 0 に収束するとか, いろいろな可能性 がある中で, 無限大に発散するというのは最悪の事が起きたという気分にもなる. このようなことが絶対に起きない方法, というのは現時点では困難な問なので, ここでの関心は 「なぜ今の方法でこのようなことが起きるのか, 特に, ある地点から狂ったように発散していく仕組 みは何なのか」をよりよく理解することにある.考察課題
前問で行った近似方法では, a = 0, b > 0 の場合に, (残念ながら) 必ず発散する解が得られるこ とを示せ.重要なメッセージは, 近似計算を気軽にやれば正しい解が得られるとは限らず, 収束・発散といっ た重要な基本的性質すら壊れているかもしれないということ, そしてそれがどのような範囲で信用 できるかは, 闇雲に計算機を使っても分からず, 数学 (両輪) が必要ということである. ヒント: 微分方程式の計算方法と, 漸化式の計算方法は酷似している事に気づいたことだろう. と いうことは, 微分方程式を解いているつもりで, その実 (少し変形された), 漸化式の計算をやってい るに過ぎないのでは. . .
10
PPM
ファイル
2 変数関数 f (x, y) を指定されたとき, その関数の値を色で表示することを目標にしよう. 1. まず講義 HP で「課題 4 で用いる題材 (間違ったプログラム集, PGM/PPM 画像サンプル)」の リンク先にある画像をダウンロードしてみよ. 2. これを画像として表示するには, display コマンドを使う. display ファイル名 3. それらのファイルをテキストエディタ mi やコマンドライン端末から less コマンドで表示して みよ. なんと, 画像でありながら中身は単なるテキストファイルであることが分かるだろう. 4. この画像ファイル形式は, グレースケールの方は, PGM (Portable Grey Map) 形式, カラーの方は PPM (Portable Pixel Map) 形式と呼ばれている. これらはプログラムで生成するのに都 合の良い, 単純な画像形式である. 5. PGM (グレー画像) の形式は以下の通り (mi で表示された結果と照らし合わせて見よ). P2 w h d G0,0 G1,0 G2,0 . . . Gw−1,0 G0,1 G1,1 G2,1 . . . Gw−1,1 . . . G1,h−1 G2,h−1 . . . Gw−1,h−1
w, h はそれぞれ横, 縦の画素数で, 整数である. d は, 最大輝度の画素–つまり白—を表す数字 で, 本演習では 255 としておく. 各 Gi,jは画素 (i, j) の値で, 0 以上 d 以下の値を指定する. 0 が 黒, 255 が白, 中間はグレーになる. 画素 (i, j) とは, 左上の画素を画素 (0,0) として, そこから 右に i 個, 下に j 個行った画素のこと. 言い換えれば, ファイル中に並んだ画素は, ちょうど横 書きの文章を書く時の字の順番で並んでいる. 6. PPM (カラー画像) の形式は以下の通り. P3 w h d R0,0 G0,0 B0,0 R1,0 G1,0 B1,0 R2,0 G2,0 B2,0 . . . Rw−1,0 Gw−1,0 Bw−1,0 R0,1 G0,1 B0,1 R1,1 G1,1 B1,1 R2,1 G2,1 B2,1 . . . Rw−1,1 Gw−1,1 Bw−1,1 . . . R1,h−1 G1,h−1 B1,h−1 R2,h−1 G2,h−1 B2,h−1 . . . Rw−1,h−1 Gw−1,h−1 Bw−1,h−1 違いは 1 行目が P3 になったことと, 各画素値を RGB の 3 つの要素で表すようにしたことだけである.
基本課題
テキストエディタ (mi など) を用いて, 10× 10 画素の, 真っ赤の PPM 画像を作り, display コ マンドで表示してみよ (表示された窓が小さすぎて, 探すのに苦労するかもしれない). ヒント: エディタでたくさん同じ行 (100 画素分) を書く事になる. こういうときはうまくコピー ペースト機能を使って, 作業を効率化すると良い. これはさすがに退屈な作業で, 次の課題のための確認作業に過ぎない. 次の課題が難なくできる 人は飛ばしてよい.基本課題
上の説明と, ダウンロードした PPM ファイルの中身をエディタで眺めて, その形式を理解せ よ. 理解したら, 512× 256 画素 (横が 512, 縦が 256) で真っ青の PPM 画像を作るプログラムを書 いてみよ. 上記形式にしたがって, データを出力する python プログラムを書き, コマンドライン端 末のリダイレクト機能を使ってデータを保存する. それを display コマンドで確認する. プログラムの概形はきっとこんな感じになるだろう. for 文の中にまた for 文があるところに注意.def print blue ppm(): print ”P3” print 512,256 print 255 for . . . for . . . print . . . print blue ppm()
11
2
変数関数の値を可視化する
PPM 形式が理解できたら, 2 変数関数を色で可視化できる.基本課題
xy 平面内の矩形領域 [−2.0, 2.0] × [−1.0, 1.0] で, 関数 f (x, y) = x 2 4 + y 2 の値を色で表示する PPM ファイルを出力するプログラムを作れ. できた PPM ファイルも提出物に含めること (ファイル名 ex11.ppm) 512× 256 の画像で上記の長方形領域全体を表すようにする. f(x, y) の値が 0 だったら黒 (0 0 0), 2.0 だったら赤 (255 0 0) ということにして, 途中の値は値に比例して補間する. 注意: 画素の値はあくまで整数である. 小数点付きの数で表示しないように注意. 小数点付きの数 x を整数にする (小数点以下を切り捨てる) のに, int(x) という関数がある. 画像を見て, 自分が縦と横を間違えていないかなどをチェックすること.12
フラクタル図形
(
マンデルブロ集合
)
関数の値を可視化できることはわかったが, 上記の関数では見ていてあまりにも退屈である. もう 少しオモローなものを表示してみよう. マンデルブロ集合は, 以下で定義される数列が n → ∞ で発散しないような複素数 c の集合で ある. z0 = c zn = z2n−1+ c (n ≥ 1) 数列 znの値を計算していく事は簡単だが, それが発散するかどうかを判定することは簡単では ない. したがってある複素数 c がマンデルブロ集合に属するかどうかを正確に判定することは難し い. ここではある打ち切り回数 N を決め, • z0, z1, . . . , zN−1の大きさ (絶対値) がどこかで 2 を超えたら発散• そうでなければ発散しない とみなすことにする (前者は一般に正しいが後者は正しいとは限らない).
オプション課題
実数 a, b, d が与えられたら, マンデルブロ集合 M を, 複素平面内の領域 [a, a + d]× [b, b + d] ((a, b) を左下隅とする一辺 d の正方形) の範囲で描画するプログラムを書け (前問の要領で PPM ファイルを出力せよ). 上で述べた方法で, ∈ M と判定された点を黒で, そうでない点を適当な色 (例えば緑) で表示せよ. 画素数は任意でよい (大きくした方がより綺麗な絵になる). 最初は a = b =−2.0, d = 4.0 とし て全体を表示し, 後は境界線付近を狙って拡大してみるとよい. おもしろい模様が見られるはずで ある. さらに余力がある人は黒と緑の 2 色だけでなく, 集合外の点をいろいろな色で表示してみよ. そ れは, z0, z1, . . . , zN−1の大きさ (絶対値) がどこかで 2 を超えた場合, 最初に超えたときの n の値に応 じて色をつけるとよい. 作品ができたら, PPM ファイルを提出物に含めてください (ファイル名 ex12-1.ppm, ex12-2.ppm, . . . )13
ビュフォンの針
基本課題
平面上に間隔 1 で無数の平行線が引かれている. その上に長さ 1 の針を適当に落とすと, 針は どのくらいの確率で平行線のひとつと交わることになるか, その近似値を計算機で求めよ. 求め方は何通りか考えられるが汎用性の高い方法は, とにかく「針を投げてみる」という試行を多数回, ひたすら計算機の中で「実験」してみる. そして何回中何回交わったかを数えることであ る. 「実験」の結果どこに針が落ちるかは乱数を用いて決定する. そのために random() とういう関 数 ([0, 1] の数を一つランダムに選んで返す) を使う. ヒント: 乱数で何を決めるかは何通りかの考え方があるが例えば以下のように考えればよい. 針の 端点を A, B とする. まず, A が最寄りの直線からどのくらいの距離に落ちるか (y) が, [−1, 1] の区 間に一様に分布する. そして, ベクトル AB が, 水平に引かれた線に対してどのくらいの角度をなす か (θ) が, [0, 2π] に一様に分布する. そこで y と θ を一様乱数で決めればよいだろう. 背景 計算機を使って, 「確率」や「期待値」を求めるのに, 乱数を使った計算は非常に汎用的な方 法論を与える. それはこの例で見たように, とにかく「試行」を乱数を用いて多数回生成し, ある事 象が起きる確率や, 試行に付随する値の期待値を求めていくという方法である. 答えを閉じた数式 として求めることができなない場合ても, 計算機を使えば数値的に答えの近似値を求めることがで きる もう少し見方を一般的にして, この方法は, どんな場合に何を計算できる手法なのかを整理して みよう. この問題で, 一回の試行は y と θ という二つの数を与えることで定まる. つまり一回の試行 は, 標本空間 Ω = [−1, 1] × [0, 2π] からの一標本とみなすことができる. Ω の部分集合として, 「針が 線と交わる (β, θ) の集合 A」が存在する. 今求めたいのは, A の面積 Ωの面積 ということである. それを求めるのに, Ω 内で点を多数発生させ, いくつが A に含まれるかでそれを計算するという のがこの手法である. もちろん, 面積と書いたのはたまたま変数が 2 つだったからで, 手法自体は変 数が 3 つでも 4 つでも通用する. その場合計算しているものは, Ω 内で定義され, fA(x) = { 1 (x∈ A) 0 (x6∈ A) で定義される関数を使って, ∫ ΩfA(x)dx ∫ Ωdx のことである. 言い換えれば, x が Ω 内で一様に発生する際の, fA(x) の平均 (期待値) を求めている. そして何も対象を fA(x) に限らなくても, この方法で任意の関数の期待値を求めることができる. ま とめると以下のようになる. 入力: 標本空間 Ω と, 関数 g(x) (x∈ Ω) 出力: x が Ω 上一様に分布する時の, g(x) の平均 ∫ Ω∫g(x)dx Ωdx
方法: 一様乱数を用いて x を Ω から多数 (N 個) 発生させ (x1, . . . , xN), 平均 1 N N ∑ i=1 g(xi) をその近似値とする. そう考えて以下の問題を解くと, 本問とまったくその原理が同じであるということに気づくこと だろう. もっともこう考えなくても, 以下は簡単に同じやり方を思いつく, わかりやすい例題だとは 思うが.
