7．栓をし忘れ浴槽に水(ルンゲ・クッタ法) 数年前のある中学校の入試

環境シミュレーション
７．栓をし忘れ浴槽に水(ルンゲ･クッタ法)
平方根に比例し，高さ h から落下する物体の落下速
｢いまさら流体力学？｣木田重雄著(丸善株式会社)より
度 2 gh と等しくなる．
数年前のある中学校の入試問題である．下図に示
す浴槽は蛇口から一定の割合で水を入れると 2 時間
で満杯になる．また，排水口は満杯の浴槽を空にす
るのに 3 時間かかる．ある日，浴槽に水を張ろうと
して，うっかり排水栓を閉じるのを忘れてしまった．
浴槽は蛇口を開いてから何時間で一杯になるだろ
う？
では，栓をし忘れて浴槽に水を張るという問題を
考えてみよう．先ず，浴槽の容積 V は 1 m 3 (底面積
S = 1 m 2 ，深さ h = 1 m)とする．問題の蛇口からは，
1 時間に浴槽全体の 2 分の 1 に相当する水が入るこ
とから，単位時間あたりの給水量(給水流量) Qin は，
浴槽容積を給水の所用時間で割った 1 2 m 3 h とな
る．この給水流量は一定である．
次に，満杯の浴槽を排水栓を抜いて空にする場合
を考えてみる．ある微小な時間 Δt の間に流れ出す
水の量は，｢トリチェリの定理｣から排出流量を
− α h とおけば，両者の積 − α h Δt となる(ここで h
は，その時点(時間 t )における浴槽の水深である)．
注意：理論的に言えば，流速が 2 gh となることから，こ
れに排出口の断面積をかけてやれば排出流量が求まるは
ずであるが，現実には，縮脈といって，流れ出す水の断面
積が排出口の断面積より小さくなる現象が現れる．また，
その小さくなる比率(縮脈率)も排水口の形状や粗度によ
って一様でないことから，理論式はそのままでは使うこと
ができない．ただ，排出流量が水深の平方根に比例すると
いう事実だけは残る．
正解？
一方，これだけの水が流れ出した結果として減少
した浴槽の水深を Δh とおけば，減少した浴槽容積
は， Δh に底面積をかけて SΔh となる．この減少し
た浴槽容積と排出された水の量は等しいはずであ
る か ら ， SΔh = −α h Δt ． 式 を 変 形 す れ ば
Δh Δt = − α h S ．ここで Δt の極限 Δt → 0 を考えれ
この問題の出題者はおそらく次のような解答を
受験生に期待していたのであろう．
この蛇口からは 1 時間に浴槽全体の 2 分の 1 に相
当する水が入る．一方，排水口からは 1 時間につき
浴槽の 3 分の 1 に相当する水がでていく．したがっ
て，蛇口と排水栓の両方を開くと 1 時間に浴槽の
1 2 − 1 3 = 1 6 に相当する水がたまっていく．したが
ば，次の微分方程式が得られる．(ちなみに微分方程
式は上記のように，浴槽への水の出入り，すなわち水とい
う物質のある系(浴槽)への収支から導き出されることか
ら｢物質収支式｣とも呼ばれる．)
って，浴槽は 6 時間で満杯になる．答え 6 時間．
一見これで良さそうに思われるが，この問題はそ
う単純ではなかった．実はこの場合｢浴槽はいつま
．．．．
でたっても一杯にならない｣が正解である．
dh
α
=−
dt
S
正解
h
ここで上式は，変数分離法で解析的に解くことが
できる．
この問題を解く鍵は，排水口から流れ出す水の速
さにある．君たちもお風呂(浴槽)の水を抜いたこと
があるだろうが，排水口から流れ出す水の速さは決
して一定ではない．浴槽が満杯に近いときは，流れ
出す水の速度は速く，空っぽに近くなるとそれにつ
れて速度も遅くなっていく．なぜそうなるかの詳し
い解説はしないが，結論だけをいうと，排出口にお
ける水の流速 v は次式のようになる．
dh
α
=−
dt
S
h → -
dh
h
=
α
S
dt → -∫ h −1 2 dh =
α
S
∫ dt
満杯の浴槽をからにするのに排水栓を抜いてか
ら 3 時間かかるという題意より，次のように α が求
まる．
v = 2 gh
0
α
H
S
− ∫ h −1 2 dh =
3
[ ]
∫ dt → − 2 h
0
0
H =1
=
α
S
[t ]
t =3
0
→ 2 1 = 3 α 1 → α = 2 3 = 0.6666... [m 5 2 h]
上式は｢トリチェリの定理｣と呼ばれ，ここでの g
は重力加速度， h は浴槽の水深である．つまり，水
深 h の浴槽から流れ出す水の速度は，浴槽の水深の
最後に，排水栓を閉じずに蛇口を開いた場合の時
- 37 -
環境シミュレーション
上式をもとに Excel で，水深 h に対して，その深
さの水を溜めるのに必要な時間 t をプロットしてみ
たのが下図である．
間 t に対する浴槽の水深 h を考えてみよう．先ほど
と同様に，ある微小な時間 Δt を考える．その間に流
れ出す水の量は − α h Δt ，また蛇口から流れ込む水
の量は，給水流量が Qin ( 0.5 m 3 h )より +Qin Δt ．そ
れら両者の結果としての，浴槽の水深の変化を Δh
とおけば，変化した浴槽容積は， Δh に底面積をか
けた SΔh となり，この変化した浴槽容積は，給水さ
れた水の量から排出された水の量を引いた値は等
しいはずである：SΔh = Qin Δt − α h Δt ．式を変形し，
Δt の極限 Δt → 0 を考えれば，
dh
dh 1 2
= Qin − α h → 1
= −
h
dt
dt 2 3
S
この微分方程式をやはり変数分離法で解いてやる．
dh
→
1 2−2 3 h
= dt → ∫
dh
1 2−2 3 h
= ∫ dt
限大となる．つまりこの問題の場合は，浴槽の水は
水深 0.5625 m までは溜まるが，それ以上は時間が
無限に経過しても溜まらない．すなわち浴槽は永遠
に満杯とはならないのである．
この場合は，さらに h = x と変数を置換してやる
ことが必要になる．積分範囲は 0 ≤ h ≤ h に対して
0 ≤ x ≤ h ．ここで h = x の辺々を微分すると
この理由は簡単である．水を張りはじめてある水
深までは水の給水流量が排出流量を上回る
( dh dt = 1 2 − 2 3 h > 0 )ために，水は溜まっていく
1 dh
1 dh
= dx →
= dx → dh = 2 xdx
2 h
2 x
これを元の微分方程式に代入して，
dt =
=(
dh
1 2−2 3 h
=
図が示すとおり，浴槽に水を張っていくために必
要時間は，浴槽の水深が 0.5 m を過ぎたあたりから
急 激 に 長 く な る ． 理 論 的 に は ，
dh dt = 1 2 − 2 3 h = 0 となる h = 9 / 16 = 0.5625 で無
2 xdx
12 xdx 9 − 3(3 − 4 x )
dx
=
=
1 2 − 2 3 x 3 − 4x
3 − 4x
9
− 3)dx
3 − 4x
が，水深が高くなるにつれて，給水流量が一定なの
に対して排出流量が徐々に大きくなっていき，ある
水深( h = 0.5625 m )になった時点で，排出流量が給水
流量と等しくなり( dh dt = 1 2 − 2 3 h = 0 )，それ以
上の水が浴槽にたまらなくなるのである．
ルンゲ・クッタ法によるシミュレーション
では同じ問題（微分方程式 dh dt = 1 2 − 2 3 h ）
辺々をそれぞれ 0 ≤ t ≤ t と 0 ≤ x ≤ h の範囲で積分
すると，以下のようになる．
h
t
を数値計算で解いてみる．ここでは差分法のひとつ
であるルンゲ・クッタ法で解いてみよう．前の章で
の少し触れたように，差分法としてはルンゲ・クッ
タ法のほうがオイラー法よりも精度がよく，実際の
コンピュータ・シミュレーションで利用されている
差分法のほとんどがルンゲ・クッタ法である．
9
∫ dt = ∫ ( 3 − 4 x − 3)dx
0
0
h
⎡ 9
⎤
ln(3 − 4 x) − 3 x ⎥
0 = ⎢−
⎣ 4
⎦0
9
9
= − ln(3 − 4 h ) − 3 h + ln(3)
4
4
9
3
∴t = ln
−3 h
4 3− 4 h
[t ]
t
ただし解説としては，サンプルプログラム全般の
説明からはじめ，ルンゲ・クッタ法については，そ
の後で説明する．
サンプルプログラムの概要
下図がこの問題のために作成したサンプルプロ
グラムの外見である．
残念ながら h(t ) = のように，時間 t の関数として
の浴槽の水深 h を求めることはできないが，逆に水
深 h に対する(水深 h の水を溜めるのに必要な)時間
t は求めることができる．
- 38 -
環境シミュレーション
次のような Sub プロシージャで定義されている．
Option Explicit
Dim qin, alpha As Double
'ﾙﾝｹﾞ-ｸｯﾀ法(Runge-Kutta Method)
Sub Sub_RKM()
Dim t, k1, k2, k3, k4 As Double
Dim i As Long
Dim h1, del_t, t_end As Double
UserForm1.TextBox1 = 0.5
UserForm1.TextBox2 = 0.1
UserForm1.TextBox3 = 10
UserForm1.TextBox4 = 0.666
UserForm1.TextBox5 = 0
このプログラムでは，B8:C9 にある｢シミュレー
ション・スタート｣というコマンドボタンを押すこ
とで，シミュレーションが開始する．すると経過時
間 t が A6 に，その時点での浴槽の水深 h が B6 に，
C6 には水槽の深さ 1 m から水深を引いた値が表示
され，その結果をうけて右の図の水色部分があたか
も浴槽中の水のように上下する．
UserForm1.Show
qin = UserForm1.TextBox1
del_t = UserForm1.TextBox2
t_end = UserForm1.TextBox3
alpha = UserForm1.TextBox4
h1 = UserForm1.TextBox5
t=0
Cells(6, 1) = t
Cells(6, 2) = h1
If MsgBox(Prompt:="Ready?", Buttons:=vbOKCancel) <>
_ vbOK Then Exit Sub
For i = 1 To t_end / del_t
t = t + del_t
k1 = dhdt(t, h1)
k2 = dhdt(t + del_t / 2, h1 + del_t * k1 / 2)
k3 = dhdt(t + del_t / 2, h1 + del_t * k2 / 2)
k4 = dhdt(t + del_t, h1 + del_t * k3)
h1 = h1 + del_t * (k1 + 2 * k2 + 2 * k3 + k4) / 6
Cells(6, 1) = t
Cells(6, 2) = h1
ActiveSheet.Calculate
If t > 0 And (h1 <= 0 Or h1 >= 1) Then Exit For
Next i
MsgBox ("ｼﾐｭﾚｰｼｮﾝ終了")
End Sub
では実際，スタートボタンを押してみよう．する
と上図のようなダイアログボックスが現れる．
ここでは， Qin = 給水流量 Qin [ m 3 / h] ，del_t =
差分法のタイムステップ Δt [h] ，t_end = シミュレ
ーションの終了時間[h]，alpha = 排出流量の速度係
数 α [m 5 2 h] ，h(0) = 時間 t = 0 での水深 h(0) [m]
Function dhdt(t As Double, h1 As Double) As Double
If h1 < 0 Then h1 = 0
dhdt = qin - alpha * h1 ^ (1 / 2)
End Function
を設定するようになっている．なお，図中ボックス
内のそれぞれの値はデフォールトとして与えてい
るものである．このようなダイアログボックスの作
成のしかたと，ダイアログボックスからの値を数値
計算へ取り込む方法については後述する．
また，この「Sub_RKM」がマクロとしてサンプ
ルプログラムの｢シミュレーションスタート｣ボタ
ンにリンクされている．（コマンドボタンへのマク
ロ登録はモンテカルロ法を参照）
またこのプログラムでは，ダイアログボックスで
パラメータを設定した後，シミュレーションの開始
前と終了後にそれぞれ下図のようなメッセージボ
ックスを表示させるようにしている．
ではプログラムの詳細をみていこう．
プログラムの詳細
Option Explicit
Dim qin, alpha As Double
'上記は，給水流量(qin)と排出流量の速度係数 α の 2 変数を倍精
'度で宣言したものである．このようにモジュールの先頭でデ
'ータ型を宣言すると，モジュール内のすべてのプロシージャが
'適用範囲となる．このような変数を｢モジュールレベル変数｣
'とよび，各プロシージャ内で宣言され，そのプロシージャ内だ
'けで有効なローカル変数と区別される．
VBA プログラム
では上記のようなサンプルプログラムを動かし
ているプログラムを見てみよう．このプログラムは
- 39 -
環境シミュレーション
t = t + del_t
'ﾙﾝｹﾞ-ｸｯﾀ法(Runge-Kutta Method)
'i がカウントされる度に，シミュレーション・クロックを del_t
'だけ進めている．
Sub Sub_RKM()
'この Sub...から End Sub までが｢Sub_RKM｣と名づけた Sub プ
'ロシージャである．
k1 = dhdt(t, h1)
k2 = dhdt(t + del_t / 2, h1 + del_t * k1 / 2)
k3 = dhdt(t + del_t / 2, h1 + del_t * k2 / 2)
k4 = dhdt(t + del_t, h1 + del_t * k3)
h1 = h1 + del_t * (k1 + 2 * k2 + 2 * k3 + k4) / 6
Dim t, k1, k2, k3, k4 As Double
Dim i As Long
Dim h1, del_t, t_end As Double
'上記 3 行は｢Sub_RKM｣内でのデータ型宣言部分である．ちな
'みに t はシミュレーション時間，k1，k2，k3，k4 はルンゲ・
'クッタ法で使用する変数，i は For...Next ステートメントのた
'めのインデックス，h1 は浴槽の水深，del_t がタイムステップ，
' t_end がシミュレーションの終了時間を表している．
UserForm1.TextBox1 = 0.5
UserForm1.TextBox2 = 0.1
UserForm1.TextBox3 = 10
UserForm1.TextBox4 = 0.666
UserForm1.TextBox5 = 0
'パラメータを設定するユーザフォーム
'にデフォールト値を入れるための行である．｢UserForm1｣とは
'そのダイアログボックスがその名称のユーザフォーム (後述
'する)に定義されていることを示す．また｢TextBox(1,2,
'3,...)｣は，数値を代入する場所がダイアログボックスで何番
'めのテキストボックス(後述する)になるかを示している．
'まず変数 k1 では，シミュレーション時間 t とその時点での水深
'h1 の値を下で定義している Function プロシージャ｢dhdt()｣
'に手渡し，その時点の dh/dt の値を返り値として k1 に受け取る．
'変数 k2 では，シミュレーション時間 t に del_t/2 を足した値
'と，さきほどの水深 h1 に del_t*k1/2 を足した値を｢dhdt()｣
'に手渡し，その時の dh/dt を返り値として k2 に受け取る．
'同様に変数 k3 では，シミュレーション時間 t に del_t/2 を足
'した値と，水深 h1 に del_t*k2/2 を足した値を｢dhdt()｣に手渡
'し，その時の dh/dt を返り値として k3 に受け取る．
'変数 k4 では，シミュレーション時間 t に del_t を足した
'値と，水深 h1 に del_t*k3 を足した値を｢dhdt()｣に手渡し，そ
'の時の dh/dt を返り値として k4 に受け取る．
'最後にルンゲ・クッタ法の仕上げとして，上記のようにして得
'られた変数 k1，k2，k3，k4 から，その平均を 1:2:2:1 の重み
'で求め，重み平均に del_t をかけた値を，そのシミュレーショ
'ンステップ内での h1 の変化量として，前ステップの h1 に加え，
'その結果を新たな h1 とする．
Cells(6, 1) = t
Cells(6, 2) = h1
UserForm1.Show
'｢UserForm1｣に定義されたダイアログボックスを画面に表示
'せよと命令している．
'更新された t と h1 を，シートのセル A6 と B6 にそれぞれ表示
'する．
qin = UserForm1.TextBox1
del_t = UserForm1.TextBox2
t_end = UserForm1.TextBox3
alpha = UserForm1.TextBox4
h1 = UserForm1.TextBox5
ActiveSheet.Calculate
'ワークシート上のグラフを描き直す．
If t > 0 And (h1 <= 0 Or h1 >= 1) Then Exit For
'ダイアログボックスに打ち込まれた数値をプログラム内のパ
'ラメータとして取り込む．
'もし(If)t が 0 より大きく(t>0)，かつ h1 が 0 以下(h1<=0)(浴
'槽から水がなくなる)かあるいは(Or)h1 が 1 以上(h1>=1)(浴槽
'から水があふれるよう)であれば(Then)，For...Next ステート
'メントを抜ける(Exit For)．
t=0
'シミュレーション・クロックを 0 に初期化する．
Next i
Cells(6, 1) = t
Cells(6, 2) = h1
'次の i で For...Next ステートメントを繰り返す．
MsgBox ("ｼﾐｭﾚｰｼｮﾝ終了")
'現在アクティブなシートのセル(6 行 1 列目つまり A6)に変数 t
'の値を，セル(6 行 2 列目つまり B6)に h1 の値を表示する．
'｢ｼﾐｭﾚｰｼｮﾝ終了｣のメッセージボックスを画面に表示させる．
If MsgBox(Prompt:="Ready?", Buttons:=vbOKCancel) <>
_ vbOK Then Exit Sub
End Sub
'｢Sub_RKM｣という Sub プロシージャの終わり．
'MsgBox(メッセージボックス)関数を続く()の条件で画面表示
'し，もし(If)それに対する回答が O.K.ボタン(vbOK)以外(<>)
'なら (Then)，この Sub プロシージャを終了する．ちなみに
' MsgBox 関数の表示条件は Prompt 文が”Ready?”で，ボタンが
'｢O.K.｣と｢Cancel｣である(先に示したメッセージボックスを見
'よ)．
Function dhdt(t As Double, h1 As Double) As Double
'｢dhdt｣という Function プロシージャの定義開始
If h1 < 0 Then h1 = 0
'さて以下の For...Next ステートメントの中がこのプログラム
'の核心部，ルンゲ･クッタ法による微分方程式の数値計算であ
'る．
'もし(If)h1 が 0 より小さい(h1<0)場合(Then)は，h1 を 0 と置
'く．
For i = 1 To t_end / del_t
dhdt = qin - alpha * h1 ^ (1 / 2)
'問題の微分方程式(dh/dt)の値を与えられた条件で計算して返
'す．
'以下の Next 行までの計算を 1 から，シミュレーション時間
't_end をタイムステップ del_t で割った回数だけ繰り返す．
- 40 -
環境シミュレーション
End Function
選択する．カーソルが｢+｣字に変わるので，フ
ォーム上でこれをドラッグして，コントロール
の配置場所とサイズを確定する．
ユーザフォームの作成方法
サンプルプログラムではシミュレーションをス
タートさせると，下図のようなダイアログボックス
（ユーザフォーム）が現れ，ここでさまざまなシミ
ュレーション条件を設定変更できるようになって
いた．
コントロールボタンの内容
[ツールボックス]にはさまざまなタイプのコントロ
ールボタンが用意されている．しかし，本授業で使
用する以下の 3 つのボタンのみである．ここでは，
それらのボタンのみについて紹介し，それ以外のボ
タンについては説明を省略する．興味のある者は，
オンラインヘルプなどで調べてほしい．
以下，Excel VBA におけるダイアログボックス
（ユーザフォーム）の作成方法を簡単に説明してお
く．
ユーザフォームの作成には先ず，[ユーザフォー
ム]シートが必要である．新規にこのシートを挿入す
るには[標準モジュール]と同様に，[Visual Basic
Editor]の[挿入]から[ユーザフォーム]を選ぶ．する
と下図のような新規の[ユーザフォーム]シートが現
れる．プロシージャと異なり，一つの[ユーザフォー
ム]シートに定義できるダイアログボックスは一つ
きりである．
１． [ラベル]：フォーム上に文字列を書き込むため
のコントロール．
２． [テキストボックス]：文字列，数値などを入力
するためのコントロール．ちなみにここに入っ
た値を VBA では，｢UserForm1.TextBox1｣の
ように参照する．
３． [コマンドボタン]：マクロ実行のボタンコント
ロール．
[テキストボッ
クス]で定義し
た入力ボックス
[ラベル]で書き込んだ文字列
上図のように[ユーザフォーム]シートで定義でき
るダイアログボックスはユーザが自由に作成・編集
できることからユーザフォームと呼ばれ，VBA の
定義済みのダイアログボックス(例えば，MsgBox
関数によるメッセージボックスなど)と区別される．
ちなみに 1 変数だけの入力を求めるのなら
「InputBox」コマンドが簡単で便利である．
この画面では，ダイアログボックスのタイトルや
サイズを変更したり，希望のボタン(正確にはコント
ロールと呼ぶ)をボックス上に配置したりすること
ができる．
１．タイトルの変更：ボックスのタイトルバーの見
出しを変更するためには，フォーム上部の文字
列(初期状態では｢UserForm1｣となっている)
を選択して，新しい名前に書き換える．
２．サイズの変更：フォームの端をクリックすると，
｢Ⅱ｣印のハンドルと呼ばれるものが四隅と四
辺上に現れるので，これをフォームが適当な大
きさになるまでドラッグする．
３．コントロールの配置：[ツールボックス]からフ
ォーム上に配置したいボタン(コントロール)を
コマンドボタンの設定方法
1) [ツールボックス]から[コマンドボタン]を選択．
「＋」カーソルをドラッグしてダイアログボック
ス上に適当な大きさのコマンドボタンを作成す
る．
2) ボタン上の文字列「CommandButton1」をクリ
ックして，これを「OK」と書き換える．書き換
えたら，ダイアログボックスの他の場所をクリッ
クして編集から抜ける．
3) 「OK」と書き換えたボタンを右クリックして，
ドロップダウンメニューから[コードの表示]を選
択する．すると次のような Private Sub プロシー
ジャが現われるので，間に「UserForm1.Hide」
を書き加え，「UserForm1」のアイコンをクリ
ックして，元のダイアログボックスのデザイン画
面に戻る．
Private Sub CommandButton1_Click()
UserForm1.Hide
End Sub
- 41 -
環境シミュレーション
ここで「UserForm1.Hide」は，同ボタンがクリ
ックされたときに，UserForm1 を隠せ（Hide）と
の意味である．
k 1 = f (t j , Y j ), k 2 = f (t j +
Δt
Δt
, Yj +
k 2 ), k 4 = f (t j + Δt , Y j + Δtk 3 )
2
2
1
1
{Y j +1 − Y j } = ( k 1 + 2 k 2 + 2 k 3 + k 4 )
6
Δt
k 3 = f (t j +
シミュレーション演習
出来上がったプログラムによって，以下のような
数値実験を行ってみよう．
ではなぜこのような差分方法(ルンゲ･クッタ法)
がオイラー法よりすぐれた O (( Δ t ) 4 ) の(累積打切り)
誤差をもつのであろうか？
１．給水流量 Qin [ m 3 / h] の変更：ダイアログボック
スで，排出流量速度係数 alpha = 0(排水口に栓
をした状態)，時間 t = 0 での水深 h(0) = 0 とお
いて，給水流量 Qin にいろいろな値を入れ，給
水流量によって浴槽が満杯になる時間がどの
ように変わるかを見てみよう．
２．排出流量速度係数 α [m
52
Δt
Δt
, Yj +
k1 )
2
2
河川の自浄モデル dL / dt = − K L をもとに，ルン
ゲ・クッタ法の局所打切り誤差が O (( Δ t ) 5 ) (すなわち
累積打切り誤差が O (( Δ t ) 4 ) )になることを示そう．
まず，自浄モデルの解析解は L ( t ) = L ( 0 ) e − Kt であ
るから，時間 t = t + Δ t におけるテイラー展開は以下
のようになる．
h] の変更：ダイアロ
グボックスで，給水流量 Qin = 0(給水がない状
態)，時間 t = 0 での水深 h(0) = 1(浴槽が満杯の
状態)とおいて，排出流量速度係数 alpha にい
ろいろな値を入れ，速度係数によって浴槽がか
らになる時間がどのように変わるかを見てみ
よう．
L (t + Δt ) = L (t ) +
+
３．栓をし忘れて浴槽に水を張りはじめた場合：デ
フォールト値のままでシミュレーションを行
い，本当に浴槽がいつまでたっても一杯になら
ないのか，浴槽の最高水位は理論値どおりかを
確認しよう．
Δ t ( 1)
( Δt ) 2 (2)
( Δ t ) 3 ( 3)
L (t ) +
L (t ) +
L (t )
1!
2!
3!
( Δt ) 4 (4)
L (t ) ⋅ ⋅ ⋅
4!
ここで，
L(1) (t ) = − KL(0)e − Kt = − KL(t )
L( 2 ) (t ) = − KL(1) (t ) = K 2 L(t )
L( 3) (t ) = K 2 L(1) (t ) = − K 3 L(t )
４．タイムステップ Δt [h] の変更：タイムステップ
del_t をいろいろ変えてみて，シミュレーショ
ンの実行速度への影響を見てみよう．
L( 4 ) (t ) = − K 3 L(1) (t ) = K 4 L(t )
より
Δt
(Δt ) 2
(K 2 )
(− K ) +
1!
2!
(Δt ) 4
(Δt ) 3
+
( K 4 ) ⋅ ⋅⋅)
(− K 3 ) +
3!
4!
L ( t + Δ t ) = L ( t )(1 +
参考書
｢理工系の基礎数学－数値計算｣高橋大輔著(岩波書
店)
次にルンゲ・クッタ法のアルゴリズムに従って演
算を行ってみる．
ルンゲ・クッタ法
ここでルンゲ・クッタ法について説明する．ルン
ゲ・クッタ法による差分方法は，オイラー法による
差分方法を
1
{Y j +1 − Y j } = f (t j , Y j )
Δt
とおけば(ここで Y j はステップ j 時点での差分解；
Y j +1 は次のステップ j+1 時点での差分解――すなわ
ち求めたい解，f (t j , Y j ) は時間 t j での与えられた微
分方程式の微分係数)，以下のように表すことができ
る．
- 42 -
環境シミュレーション
f (t , L(t)) = − KL(t )
k1 = f (t , L(t)) = − KL(t )
Δt
Δt
Δt
, L(t) +
k1 ) = − K ( L(t) +
k1 )
2
2
2
Δt
Δt
(− KL(t ))
k1 = − KL(t) − K
= − KL(t) − K
2
2
Δt
= (− K + K 2 ) L(t )
2
Δt
Δt
Δt
k3 = f (t + , L(t ) +
k2 ) = − K ( L(t ) +
k2 )
2
2
2
Δt
Δt
Δt
k2 = − KL(t ) − K
(− K + K 2 ) L(t )
= − K L(t ) − K
2
2
2
2
t
(
t
)
Δ
Δ
) L(t )
= (− K + K 2
− K3
2
4
k4 = f (t + Δt , L(t) + Δtk3) = − K ( L(t) + Δtk3 )
k2 = f (t +
= − KL(t) − KΔtk3
(Δt )2
Δt
) L(t )
− K3
2
4
(Δt )2
(Δt )3
) L(t )
= (− K + K 2 Δt − K 3
+ K4
2
4
Δt
Δt
k1 = − K
L( t )
6
6
Δt
( Δt ) 2
2 Δt
) L(t )
k2 = ( − K
+ K2
6
3
6
2 Δt
Δt
( Δt ) 2
( Δt ) 3
+ K2
− K3
) L(t )
k3 = ( − K
6
3
6
12
= − KL(t) − KΔt (− K + K 2
( Δt ) 2
( Δt ) 3
( Δt ) 4
Δt
Δt
+ K2
− K3
+ K4
) L( t )
k4 = ( − K
6
6
6
12
24
Δt
( k1 + k2 + k3 + k4 )
L( t ) +
6
= L ( t )(1 +
+
(Δt ) 2
(Δt ) 3
Δt
(− K ) +
(K 2 ) +
(− K 3 )
1!
2!
3!
(Δt ) 4
( K 4 ))
4!
となり，確かにルンゲ・クッタ法の局所打切り誤差
は O (( Δ t ) 5 ) となる．
別の角度からルンゲ・クッタ法の局所打切り誤差
が O (( Δ t ) 5 ) となることを示そう．下の表は，5 つの
関数 f ( t ) = t , t 2 , t 3 , t 4 , t 5 について， f (1) の値を，初
期値 t = 0 , Δ t = 1 で，3 つの差分(オイラー，ホイン，
ルンゲ・クッタ)法で求めてみたものである．
表よりも，ルンゲ・クッタ法では f ( t ) = t ～ t 4 ま
で誤差なく計算できるが， f ( t ) = t 5 より高次な関
数ではエラーが生じることがわかる．すなわちこの
ことこそ，ルンゲ・クッタ法の局所打切り誤差が
O (( Δ t ) 5 ) となることと同値である．
- 43 -
環境シミュレーション
オイラー，ホイン，ルンゲ・クッタの比較
f (t )
t
t2
t3
t4
t5
f ′(t )
1
2t
3t 2
4t 3
5t 4
0 + 1 ⋅ (1) = 1
0 + 1 ⋅ (2 ⋅ 0 ) = 0
0 + 1 ⋅ (3 ⋅ 0 2 ) = 0
0 + 1 ⋅ (4 ⋅ 0 3 ) = 0
0 + 1 ⋅ (5 ⋅ 0 4 ) = 0
オイラー
ホイン
0 + 1⋅
1
(1 + 1) = 1
2
ルンゲ・クッタ
0 + 1⋅
1
{
6
1+1+
4 ⋅ (1)}
=1
0 + 1⋅
0 + 1⋅
0 + 1⋅
0 + 1⋅
1
( 2 ⋅ 0 + 2 ⋅ 1) = 1
2
1
3
( 3 ⋅ 0 + 3 ⋅ 1) =
2
2
1
4
( 4 ⋅ 0 + 4 ⋅ 1) =
2
2
1
5
( 4 ⋅ 0 + 5 ⋅ 1) =
2
2
0 + 1⋅
1
{
6
2 ⋅ 0 + 2 ⋅1 +
0 + 1⋅
1
{
6
3 ⋅ 0 2 + 3 ⋅1 2 +
0 + 1⋅
1
{
6
4 ⋅ 0 3 + 4 ⋅13 +
1
{
6
4
5 ⋅ 0 + 5 ⋅1 4 +
1
4 ⋅ ( 2 ⋅ ( ) )}
2
=1
1
4 ⋅ ( 3 ⋅ ( ) 2 )}
2
=1
1
4 ⋅ ( 4 ⋅ ( ) 3 )}
2
=1
1
4 ⋅ ( 5 ⋅ ( ) 4 )}
2
100
=
96
0 + 1⋅
最後に，Sub プロシージャを用いず，全てをワークシート上でシミュレーションするプログラム例を示し
ておく．
G3
=G2+$B$5/6*(H2+2*I2+2*J2+K2)
H3
=$B$1-$B$2*G3^(1/2)
I3
=$B$1-$B$2*(G3+$B$5*H3/2)^(1/2)
J3
=$B$1-$B$2*(G3+$B$5*I3/2)^(1/2)
K3
=$B$1-$B$2*(G3+$B$5*J3)^(1/2)
- 44 -