授業資料/13 の変更点

追加された行はこの色です。
削除された行はこの色です。
授業資料/13 へ行く。
授業資料/13 の差分を削除
#contents

* 前回の課題について [#z81ce357]
前回の総仕上げは，与えられた地図の二点間で最短経路を道が何通りあるか数えよというものであった．
(地図については前回の web を参照のこと)

図が非常に素直で簡単な形をしているので，最短経路の個数を数えるには単純なルールでできることは前回のヒントからわかったと思う．
具体的には，
&br;
「ある点を通る最短経路は，その上の点か右の点を通ってやってくる」
&br;
ので，これは言い換えると，
&br;
「ある点を通る最短経路の数 =  その上の点を通る最短経路の数 + その右の点を通る最短経路の数」
&br;
ということになる．あとはこれを具体的にしていけばよい．
例えば地図を素直に x,y 座標で以下のように特徴付けたとして，
CENTER:&ref(./route-exp.png,70%);

のように考えたとしよう．
このとき，座標 (i,j) (0 &le; i &le; 6, 0 &le; j &le; 4)の点を通る最短経路の数を f(i,j)とすると，数式に直せば以下のような関係式が成り立つことが分かる．
CENTER:&ref(./map-eq-1.png);

となるので，これを使えばよい．あとは，この関係式を答えに近い方から使うか遠い方から使うかでプログラムの書き方が違うだけだ．

なお，上の式を途中まで手で展開した形の
CENTER:&ref(./map-eq-2.png);

を使っても良い．ただ，こういうように「人間が途中まで介入する時に計算を失敗することがある」ので，一般的にはどこまで介入するかは悩ましい問題だ．この問題の場合はまず間違えないだろうが…

&br;
** 再帰定義を使う場合 [#q58a3e88]
これは数式を素直に使えばよいので，例えば以下のようになるだろう．
// programu source 表記
#highlighter(language=ruby,number=on,cache=on){{
include Math

# f を求める再帰定義式
def f(i,j)
  if ((i==6) && (j==4)) then
      return 1
    elsif (j==4) then
      return f(i+1,4)
    elsif (i==6) then
      return f(6,j+1)
    else
      return f(i+1,j) + f(i, j+1)
  end
end

# 以下，プログラム本体

# f(0,0) の値を出力する．
print(f(0,0),"\n")
}}

&br;
** 動的計画法もどきを使う場合 [#ge33c2ea]
右上の Start 地点から，「何ステップ目か」ということに着目すれば，ステップ数を増やしていくように計算を進めればよいことに気づく．つまり，次の図の赤い文字の順に計算が進められるべきだということがわかるだろう．
CENTER:&ref(./route-exp-2.png,70%);

少し工夫がいるが，まあ例えば次のようにプログラムを書けるだろう．
// programu source 表記
#highlighter(language=ruby,number=on,cache=on){{
include Math

# step の小さいほうから計算していって，Goal での値を返す. 一応，地図のサイズを一般にしておこう.
def calc_top(nx,ny)

  # 各点を通る最短経路数を記録しておく配列を用意しよう．
  f = Array.new(nx+1){ Array.new(ny+1){ 0 } }

  # Start 地点での f は 1 だね．
  f[nx][ny] = 1

  # 総ステップ数
  all_step = nx + ny

  # step 1 から all_step まで
  for step in 1..all_step do

    # 同じ step の点は x + y =  all_step - step を満たすことを利用して…
    for i in 0..(all_step - step) do
      j = (all_step - step) - i
      if ((0 <= i) && (i <= nx) && (0 <= j) && (j <= ny)) then # 地図に入るときだけ考慮すればよい

        # あとは普通に計算する
        if (i==nx) then
            f[i][j] = f[i][j+1]
          elsif (j==ny) then
            f[i][j] = f[i+1][j]
          else
            f[i][j] = f[i+1][j] + f[i][j+1]
        end

      end
    end

  end

  # f(0,0) を返す．
  return f[0][0]
end

# Goal での f の値を計算する
ans = calc_top(6,4)

# 知りたい値を表示する
print(ans,"\n")
}}

** 少し賢い考え方 [#ec4aa7e8]
最短経路数の依存関係を考えると「右上から左下に斜め」に話が進むことがわかる．
言い換えると，「この斜めの方向が素直な方向」になるので，この方向に図を傾けて考えるのが論理的に素直だということになる．
実際に傾けてみよう．
すると，依存関係が前回の「木構造」のところでやった話と全く同じなので，じつは同様に「木構造」で考えるとより分かりやすいということも推測がつく．
そこで，たとえば全体が「わかりやすい木構造」になるようにさらに地図を少し描き足し，全体のつじつまが合うように「葉」に数字をつけると以下のような図になる．
CENTER:&ref(./route-tree.png,60%);

あとは，下の数字を足すと自分の数字になる，という仕組みを前回と同じようにプログラムするだけだ．
これで考えると，動的計画法もどきのプログラミングも上と違って随分簡単に書けるはずだ．
計算量は増えるが，こうした素直なプログラムにはミスが入る確率は段違いに低くなる．
答えが間違っていれば意味が無いのだから，まずはこうした「素直な」「確実な」プログラムを目指そう．
&br;
このように，
&br;
CENTER:&size(24){''依存関係などをよく考えて，''};
CENTER:&size(24){''事前に問題そのものを処理しやすい形に理解し直すことが''};
CENTER:&size(24){''プログラミングには重要である．''};
&br;

* 総合力: 微分方程式を解いて, 解を見てみよう [#u0bdf45e]
これまででプログラミングについて''最低限''の知識は身につけたので，少し応用へ向けて実習を行ってみよう．
いろんな題材があるが，まずはごく簡単な常微分方程式に対して，近似計算をして近似解を求めてその様子をみるという一連の流れをやってみよう．

&br;
** 自由落下をシミュレートしてみる [#w4302938]
学生が最もよく知っている常微分方程式である，Newton 力学に基づく運動方程式を題材に取り上げよう．
まずもっとも簡単な，「真空中を物体が自由落下してくる様子」を考えよう．

時間を t sec, 地上からの物体の高さを u(t) m，重力加速度 g = 9.8 m/(s^2) とすると，
Newton の運動方程式 F = ma より，
CENTER:&ref(./newton.png,70%);
CENTER:''物体の高さについての常微分方程式''

という常微分方程式が得られる．
これは二階常微分方程式で少し扱いにくいので，速度 v = du/dt という新しい従属変数を導入して，
CENTER:&ref(./newton-2.png,70%);
CENTER:''物体の高さについての常微分方程式(変形版)''

というものに変形して，これについて考えよう．
この方程式は手で解けてしまうが，手で厳密解が求まらない場合にも通用する考え方として今回はプログラムで近似解を求めてみよう．それには様々な手法があるが，今回は最も単純な ''Euler法'' というものを使ってみよう．

&br;
*** Euler 法 [#q2524109]
Euler 法というのは簡単で，時間発展項を前進差分で近似することで，古い時間の変数値を代入すると新しい時間の変数値が出てくる関係式に問題を変形するものである．

まあ，理論を言うより見たほうが早いだろう．
今回の場合，時間変数 t を △t の幅で「刻む」ことにして，上の常微分方程式を
CENTER:&ref(./euler.png,70%);
CENTER:''上の常微分方程式を Euler 法で「近似」したもの''

とするのが Euler 法である．そしてこの式をきれいに書き直すと，
CENTER:&ref(./euler-2.png,70%);
CENTER:''上の常微分方程式を Euler 法で「近似」したもの(左辺が「新値」，右辺が「旧値」)''

となるので，これをみると「右辺に知っている値を代入するだけで，新しい時間での値が得られる」形になっている．
常微分方程式をこのように近似式に置き換えるのが Euler 法である．
これならプログラムが書けそうだという感触が得られるだろう．
&br;

&ref(/materials/warning.png); もちろん，Euler 法で作った「近似式」を解いて得られるのは「近似解」であり，もとの常微分方程式の解ではない．ただ，そんなに違わないんじゃない？ という期待をするということだね．
&br;

*** 問題設定 [#o50d4c8e]
さて，シミュレートするには具体的な条件を詰めないといけない．
そこで，以下のように初期条件などを設定しよう．

''[初期条件]''
初期高さ u(0) = 1000.0 (m)
初期速度 v(0) = 0.0 (m/s)

''[終了条件]''
高さが 0 m以下となったらシミュレーション終了.

''[近似条件]''
時間刻み幅 △t = 0.1 sec

&br;
*** プログラム [#d15c251c]

// 実習アイコン
&ref(/materials/notes.png); 以下の部分は実習しながら理解しよう．
&br;

プログラムは上の条件に従って計算して 時間, 高さ の数字ペアを出力するものとしよう．
すると，だいたい以下のような感じになるだろうか．
// programu source 表記
#highlighter(language=ruby,number=on,cache=on){{
include Math

# Euler 法に基づいて，次の時間ステップの u,v を計算する
def euler(u,v, dt)
  g = 9.8

  u_new = u + v*dt
  v_new = v - g*dt

  return u_new, v_new
end

# 時間と高さを出力する
def output(u, time)
  print(time,", ",u,"\n")
end

# 以下，プログラム本体

# 初期値
u = 1000.0
v = 0.0
time = 0.0

# 時間刻み幅
dt = 0.1

# 高さ u が正のうちは計算を繰り返す
while (u > 0) do

  # 高さと時間を出力
  output(u, time)

  # Euler 法で次の時間ステップの u,v を計算する
  u,v = euler(u,v,dt)

  # 時間を更新
  time += dt
end
}}

&br;
*** プログラムを動かしてデータを得る [#r786047b]
上のプログラムを実行して，データを得よう．得られたデータが画面に表示されるだけではあまり役に立たないので，ファイルに保存しよう．
それには例えば，上のプログラムファイル名が euler.rb だとして，
  ruby -w euler.rb > euler.dat

などとして実行すればよい…のは覚えているね？ この場合，データが euler.dat というファイルに保存される．

&br;
*** 得られたデータを動画的に見てみよう [#h76da810]
上のデータファイルを見てみると，全部で 144行ある．そこで，この 144行を 1行ずつ gnuplot などで見てみれば動画的に状況が見られることになる．
そこで，以下のように作業してみよう．
&br;

1) gnuplot の動作ファイルを作る:  ''euler.plt'' という名前のファイルを作り，次の中身を書き込む．
// programu source 表記
#highlighter(language=ruby,number=on,cache=on){{
set xrange[0:15]
set yrange[0:1000]
set xlabel "Time"
set ylabel "Height"
unset key

line=0
load "euler.loop"
}}
&br;

2) gnuplot の動作ファイルを作る(1行ずつ表示する部分): ''euler.loop'' という名前のファイルを作り，次の中身を書き込む．
// programu source 表記
#highlighter(language=ruby,number=on,cache=on){{
plot "euler.dat" every ::line::line w p pt 6 ps 3
line=line+1
pause 0.2
if(line<144) reread
}}
&br;

3a) gnuplot (メニュー → Math → gnuplot) を起動し，
3b) 上のデータファイル ''euler.dat'', gnuplot の動作ファイル ''euler.plt'', ''euler.loop'' の3つのファイルが置いてあるディレクトリに gnuplot の動作ディレクトリを「移動」させる(「移動」とか「cd」とか書いてあるボタンがあるだろう)．
3b) 上のデータファイル ''euler.dat'', gnuplot の動作ファイル ''euler.plt'', ''euler.loop'' の3つのファイルが置いてあるディレクトリに gnuplot の動作ディレクトリを「移動」させる(「移動」とか「ChDir」とか書いてあるボタンがあるだろう)．
3c) gnuplot の命令画面で
  load "euler.plt"

と入力する．問題がなければ，以下の図のような表示で動画的にその高さを示す ○ 印が動いていくはずだ(横軸に時間をとっている)．
CENTER:&ref(./gnuplot-movie.png,60%);
&br;

// 注意アイコン
&ref(/materials/warning.png); gnuplot で1行ずつ表示を繰り返させるのに使っているテクニックは主に以下の二つである．
- ''reread'' : これは，その命令が入ったファイルを読み直すという意味である．if と組み合わせて Ruby の while みたいに使っている．
- plot 命令の後ろについている ''every ::line::line'' : これはデータの第 line+1 行目のみを読んで plot しろということを意味する．
&br;

&ref(/materials/warning.png); すべての瞬間の画像を用意し，それをつなげて動画を作るという方法もある．余裕のあるものはこちらにチャレンジするのもよいだろう．ちなみに，画像をつなげて動画にするソフトウェアとしてWindows 系ならば

- ウィンドウズムービーメーカー
(大学のPC の場合: メニュー → プログラム → アクセサリ → ファイル名を指定して実行 → moviemk.exe と入力 で起動): Windows に付属の無料ソフト
- [[TMPGEnc 無料版:http://www.tmpgenc.net/ja/j_download.html]]

などが，Unix 系(含む Mac)ならば

- ffmpeg

などがあるので(自宅でインストールが必要なものはインストールして)使ってみよう．
&br;

* 今日の総仕上げ [#d01db1cd]
上の例同様に，今度は重力下で Newton の運動方程式に沿って動く「振り子の運動」をシミュレーションしてみよう．
CENTER:&ref(./pendulum.png,60%);

まずは各条件を以下のようにしてみよう．

''[状況条件]''
振り子の紐の長さ L = 1.0 m

''[初期条件]''
初期角度 &theta;(0) = &pi;/3 (rad) = 60.0 (度)
初期速度 d&theta;/dt(0) = 0.0 (rad/s)

''[終了条件]''
t = 0 sec から開始して，各自自分で決めた時間 t = T_limit までシミュレーションを行う.
t = 0 sec から開始して，各自自分で決めた時間 t = T_limit (5.0ぐらい?)までシミュレーションを行う.

''[近似条件]''
時間刻み幅 △t = 0.1 sec
時間刻み幅 △t = 0.01 sec (少し小さくしておいた)

''[近似方法]''
Euler 法
&br;

&ref(/materials/warning.png); 今度は空間が二次元の問題だが，時間発展する変数は「角度」一つだけに帰着できる．
そこで，まずは Newton の運動方程式を角度に対する常微分方程式に変形しよう．
それから，上の例のように，従属変数を新たに加えることで問題を 1階の常微分方程式に変形しよう．
そうすれば Euler 法で近似的に解けるだろう．
あと，プログラムのデータ出力時に x, y の二次元の座標に変換しないと動画が見えないので，そこは気をつけよう．
&br;

** 上級者向け(1) [#fc3e4a8a]
「二重振り子」(一つめの振り子にさらに振り子がぶら下がっている)について同様にシミュレーションしてみよう.

&br;
** 上級者向け(2) [#rbd96073]
Euler 法ではなく，古典 Runge-Kutta 法とよばれる方法について調べ，これを近似として使ってみよう．


&br;
* レポート [#v61f4984]

本日受けた講義および行った実習について，簡単にまとめた報告をせよ.
また，総仕上げの部分で自分で書いたプログラムを報告せよ.

もちろん各自の

+ 所属(学部，学科)
+ 学籍番号
+ 学年
+ 氏名
+ 日時
+ 肝心のレポート内容(得た知見，作業について気づいたこと等も)

を書くのを忘れないように. 

* about Icons, ClipArts [#d80975f6]
Some icons in this page are downloadable at [[ICONFINDER:http://www.iconfinder.net/]].

The "note" icon &ref(/materials/notes.png); designed by [[Marco Martin:http://www.notmart.org/]] is distributed with the LGPL licence,
the "warning" icon &ref(/materials/warning.png); designed by [[Alexandre Moore:http://nuovext.pwsp.net/]] with the GPL licence
and the "triangle" icon &ref(/materials/JNorth_arrow-right-sm.png); designed by [[Joseph North:http://sweetie.sublink.ca/]] is distributed with the [[Creative Commons (Attribution-Noncommercial-Share Alike 3.0 Unported):http://creativecommons.org/licenses/by-nc-sa/3.0/]] licence.

Some clip arts used in this page are downloadable at [[Open Clip Art Library:http://www.openclipart.org/]].
We deeply appreciate their superb works. With licence, they describe that "the actual clipart content on open clipart library is Public domain" in the web.

// ━┃┏┓┛┗┣┳┫┻╋


// コマンドライン入力は「行頭をブランクで始める」.
// コマンドライン出力は「行頭を > で始める」.

// 実習アイコン
// &ref(/materials/notes.png); 

// 注意アイコン
// &ref(/materials/warning.png);

// Link アイコン
// &ref(/materials/JNorth_arrow-right-sm.png);

// OK アイコン
// &ref(/materials/OK.png);

// NG アイコン
// &ref(/materials/NG.png);

// 大文字での強調 
// CENTER:&size(24){''ほげほげ''};

// programu source 表記
// #highlighter(language=ruby,number=on,cache=on){{}}