気候システム論

ふろおけモデルの簡単な例、そのRによる計算プログラム

【[2015-06-07改訂] 表計算ソフトウェアによる計算シート例を [別ページ]として追加しました。】

[別ページ]で述べた「ふろおけモデル」のうちでも簡単な場合として、流入量(単位時間あたりの質量)F_inは一定値であり、流出量(単位時間あたりの質量)F_outが水位 h に比例する場合を考える。ここでは水平断面積 A と水の密度 ρ は一定値と仮定するので、ふろおけにたまっている水の質量(以下、便宜上「ストック」と呼ぶ) S = ρ A h をおもな変数にとって、「流出量F_outはストックSに比例する」と考えることにする。式では次のように書ける。係数 c (便宜上「流出の係数」と呼ぶ)は定数である。

dS/dt = F_in - F_out ... (1)

F_out = c S ... (2)

(1), (2)をあわせてF_outをひとまず消去し、 Sは時間tの関数、F_inとcは定数であることに注意すると、このモデルは次のようなSに関する常微分方程式で表現できる。

dS/dt = F_in - c S ... (3)

ここで、仮に dS/dt = 0 という条件を置いてみると、それを満たす解(定常状態 steady state)がある。添え字「steady」をつけて示すことにする。 F_out,steady = c S_steady = F_inであることから、 S_steady = F_in / c である。これを使って、(3)を S - S_steady = S' についての式に書きかえると

dS'/dt = -c S' ... (4)

となり、その解は、指数関数 exp(x) = e^x を使って

S'(t) = S'(0) exp( -c t ) ... (5)

のように書くことができる。ただしS'(0)は時刻 t = 0 でのS'の値である。Sについての式にもどすと

S(t) = S_steady + (S(0) - S_steady) exp( -c t ) ... (6)

となる。cは正の定数なので、(6)の右辺第2項は、時間 t が大きくなると 0 に近づく(減衰する)。つまり、初期値S(0)がどんな値であっても、Sは時間とともにS_steadyに近づく。

このように、この例では、微分方程式が数式の形で解ける(「解析解が得られる」と表現されることが多い)。

時間微分の差分近似

モデルが少し複雑になると、解析解が得られなくなる。その場合に有効なのは、数値計算によって近似解を求めることである。有限の時間区間(時間差Δt)の間のSの変化 ΔS = S(t+Δt) - S(t) を考える。 dS/dtは ΔS/Δt の Δt→0 での極限だが、これを、極限をとらないΔS/Δtで近似する。たとえば (1) の近似は次のようになる。

ΔS/Δt = F_in - F_out ... (7)

これを変形して、S(t)がわかっているときにS(t+Δt)を求める式が得られる。

S(t+Δt) = S(t) + ( F_in - F_out ) Δt ... (8)

いまF_inは定数としたが、F_outは時間とともに変化する。 (8)式を使おうとすると、時間区間のうちのどの時刻でのF_outの値を使うかが問題になる。ここでは、区間の初めの時刻 t での値で代表させることにする。

【微分を有限区間の差分で近似する方法にはいろいろなものがある。気象学で応用される方法の概略は、別の教材ページ「数値計算の教材: 時間発展問題の差分計算」で紹介した。区間の中央で代表させる方法もあり、うまく使えればそのほうが精度がよいのだが、数値があばれること(「計算不安定」という)をふせぐのにくふうがいるので、ここでは、区間の初めで代表させる(専門用語で「前方差分」という)方法に限って紹介する。】

Δtが大きすぎると、差分による計算結果は微分方程式の解の近似として不適切になる。 Δtを小さくすると、計算量がふえて、計算機に対する負荷が重くなるほか、誤差の累積が問題になる場合もある。適切なΔtを選ぶには、試行錯誤を必要とすることもあるが、まず重要なのは、起こる現象の時間規模よりもやや細かくとることである。たとえば、(6)式のような指数関数型の変化が起こるならば、Δtは 1/c よりも小さく(たとえばその5分の1くらいに)する。

R言語によるプログラム例

わたしが教材として示すプログラム例は、百行くらいまでの分量のソースプログラム(人が読み書きする形でのプログラム)を必要に応じて部分的に書きかえて実行することができる人向けである。プログラミングをしない人に使ってもらうために配慮したプログラムではない。また、計算のプロフェッショナルの観点ですぐれたプログラムでもない。

わたしはこれまで、数値計算のプログラムをおもにFortran言語で書いてきた。「数値計算の教材: 時間発展問題の差分計算」の例文は、Fortran, Pascal, C言語によるものを用意した。

ここでは、R言語によるものを整備したいと考えている。 Rの処理系は、農工大のパソコン教室のWindowsマシンにはすでにはいっているし、オープンソースソフトウェアとして提供され、 MS Windows, Mac, Linuxをはじめとする多くのOSの計算機上で動くので、 Rに慣れることは多くの進路の人にとって損にならないと思う。 Rは統計的データ解析のための言語であり、ここで行なうのは数値シミュレーションの基礎となる理論的計算なので、 Rの典型的な使いかたではないのだが、なんとか使うことができる。

なお、Rを含むほとんどのプログラム言語は、数値を扱うことはできるが、単位のついた物理量を扱う機能をもっていない。物理量の次元と単位についてはプログラミングの前の段階で人が注意して扱わなければならない。

教材用のプログラムを作るにあたって、次のような方針をとることにする。

定数値、初期値、計算対象時間の長さ、表示する結果の選択などは、プログラムを書きかえる形で与える。少数箇所を書きかえればすむようにくふうしておく。 (プログラムを書きかえずに多くの種類の計算が実行できるような汎用的構成にはしない。)
今回の例ではF_inは定数だが、それが変化しうる場合へのプログラムの書きかえがしやすいように、変数として表現し、そこに定数を代入する形にする。
結果をグラフ表示するつごうで、計算結果を「配列変数」にとっておく形にした。あらかじめ、時間ステップをいくつ計算するか決めて、その寸法の配列変数を確保する。 (結果表示をかまわなければ、ひとつの変数の値を毎ステップ書きかえる形のプログラムも可能であるが。)
グラフの作図は、計算しながら結果を確認する目的だけを考え、報告のためのくふうは(まだ)しない。パッケージggplot2かlatticeを使うのがよいらしいが、わたしがこれまでに習得できたのが R処理系常備の(パッケージ使用をことわらなくても使える) plot その他の機能の初歩的な使いかたなので、今回はとりあえずそれを使った例を示す。(今後の改訂で、別のものにさしかえる可能性がある。)
モデル内の数量の単位はSI単位とする。数量を別の単位で考えた場合は、SIへの換算を機械にやらせるようにプログラムする。
数量の規模は、[ふろおけモデルのページ]で文献としてあげた MEADOWSの本の最初のふろおけの例題 (数理的構造は今回の問題と違うが)を大まかな参考にして決めた。

プログラム例を、下の「ソースプログラム」という見出しのあとに用意した。その部分だけの内容をテキストファイルとして、たとえば bathtubex1.R という名まえで保存したとする。 R処理系を起動し、 Windows用の処理系の日本語メニューで言えば「ファイル」の中の「Rコードのソースを読み込み」でそのファイルを指定すれば、実行されるだろう。しかしむしろ、「ディレクトリの変更」をしたあとで、「R console」のコマンドとして(「>」のプロンプトから) 「source "bathtubex1.R"」のような形で実行するのがよいだろう。プログラムの編集は、Rの外のテキストエディタ (Windowsの「メモ帳」を含む)ですることもできるが、メニューの「ファイル」の中の「スクリプトを開く」からRのエディタを使って「保存」するのが便利だろう。

プログラムのどこを変えたら結果がどう変わったかを確認するには、プログラムの複数のバージョンを比較できる形で持つ必要がある。エディタで単純に保存すると上書きになってしまうので、変更前のプログラムをコピーしておいてファイル名で区別するなどのくふうを、各自してほしい。 (バージョン管理を支援するソフトウェアを使うとよいのかもしれないが、わたしは(まだ)使っていない。)

Rのプログラム中のコメントや文字定数などに日本語文字を使ってもかまわないはずだが、文字コードの問題にぶつかるおそれがあるので、例文中の文字はASCIIに限り、コメントは(日本語ローマ字でもよいのだが、ここでは)英語でつけた。

R言語では各行の「#」という字からあとはコメントであり計算処理上は無視される。

例文のプログラムを部分に分けて説明するために、コメントに(1)(2)のような番号を入れた。(式の番号ではない。)

(1) モデルのパラメータ

流入量 F_in (変数名inflowc)と、流出の係数 c (変数名outcoef)を与える。 inflowc の値は 10 kg/分と考えた。 outcoef は直観的に与えにくいので、この例に限っては解析解がわかっていることを利用し、定常状態となるストック量 S_steady (変数名stockstd)のほうを 50 kgと決め、それが定常状態となる条件から計算された流出の係数の値をひとまずoutcoefstdとする。例文ではその値をそのまま outcoef に代入している。

R言語で「<-」は代入をあらわす記号である。右側の式の計算結果を左側の変数にいれる。

(2) シミュレーション実行のパラメータ

差分の1ステップの時間間隔Δt (変数名deltat)を1分、シミュレーションのステップ数(変数名 nstep)を20 (つまり継続時間を20分)とする。 (シミュレーションの考えかたのうえでは継続時間から考えて割り算でステップ数を決めるべきだが、割り算の結果が整数にならないとやっかいなので、この手順にした。) また、初期(t = 0)でのストック量S(0) (変数名 stockinit) を 30 kgとしてみる。 (これをS_steadyと同じにするとつまらない結果になるので、それと同じ桁で別の値を入れてみる。)

(3) 配列変数を確保する

時間とともに変化する変数(および、この例では一定だが、別の問題では変化させたい、流入量inflow)について、それぞれ、(nstep+1)個の記憶場所を確保する。この手続きの書きかたは、プログラム言語によって違う。なお、(3x)と注記した部分は、(5)の理論値(theoretical values)の計算のためのものなので、ふつうのシミュレーションの場合は、なくてよい。

(4) 差分法による計算

シミュレーションの本体。まずストック量の初期値を与えたあと、 for文によるループで、時間をΔtずつ進めながら、 nstep回、次の時間ステップでのストック量を計算する。

【[2015-06-09補足] この部分のプログラムは、R言語の「for(k in 1:length(time))」のループで、 kの値を1から順にひとつずつ length(time) [この例文では21] まで変えてそれぞれ実行されることを前提としている。もし並列に、あるいは順不同に実行される場合は期待どおりの動作をしない。その場合は確実に順に実行されるようにプログラムを変更する必要がある。】

(5) 理論値(解析解)の計算

この例題には、上に述べたように、解析解があるので、それも計算してみる。 (解析解があれば差分計算する必要はないのだが、両者の結果を比較することによって差分計算の方法の精度を評価することができ、解析解がない問題についての差分計算の結果を評価する際に役だつ。)

(6) 結果を表示する

ここでは、ストック量の時系列を表示する手続き(6a)と、流入量・流出量の時系列を表示する手続き(6b)を用意した。両者の切りかえは、(条件つき選択の構文などを使ったプログラムを組むのがよいと思うが、とりあえず)、ソースプログラムを編集して、行の頭にコメント印「#」を打って命令を無効にするか、それを消して命令を有効にするか、という形によっている。

数値を数字で表示する手続きも用意したほうがよいのだが、プログラム実行直後に、Rの処理系内のコマンド行として(プロンプト「>」のあとに) 「stock」や「outflow」などの変数名を単独で打てば、その変数の値が表示される (配列の場合には列挙される)ので、この例題の結果を確認することに限ったとりあえずの目的にはそれで充分かと思う。

実験の勧め

そのまま実行するほかに、(1), (2)の部分を次のように変更して結果を比較してみてほしい。

流入量「inflowc <- 10/60」の右辺を変える。たとえば、2倍(20/60)にしたり、半分(5/60)にしたり、0にしたりする。
流出量の係数「outcoef <- outcoefstd」の右辺を変える。たとえば、2倍 (2*outcoefstd)にしたり、半分 (0.5*outcoefstd)にしたりする。
ストック量の初期値「stockinit <- 30」の右辺の値を変える。
シミュレーションの時間ステップ数「nstep <- 20」の右辺の値を変える。時間間隔 deltat を変えなければ、シミュレーション対象時間の長さを変えることになる。

【[2015-06-09暫定] ひとつのプログラムの中で複数のパラメータセットによるシミュレーションを実行し、結果を比較できる形で表示できるようにしたいが、そのためのプログラムの変更のしかたの説明は、まだ用意していない。追って用意したい。】

ソースプログラム

# bathtubex1.R
# A bath-tub model
# A simple case: outflow is proportional to stock
# R program 2015-06-02, latest revision 2015-06-04 MASUDA Kooiti

##### (1) parameters of the model
inflowc    <-  10/60            # inflow rate, constant (kg/s)
                                # 10/60 ... 10 kg/min
stockstd   <-  50               # steady-state level of stock (kg)
outcoefstd <-  inflowc/stockstd # steady-state level of outflow/stock (1/s)
outcoef    <-  outcoefstd       # outflow/stock (1/s)

##### (2) parameters of a run of simulation
deltat  <-  60          # time interval of a step (s)
nstep   <-  20          # number of time steps
stockinit <-  30        # initial level of stock  (kg)

##### (3) preparing array variables
time    <- seq(0, deltat*nstep, deltat)
stock   <- time * 0  # stock (kg)
inflow  <- time * 0  # inflow rate  (kg/s)
outflow <- time * 0  # outflow rate (kg/s)
##### (3x) (array variables for theoretical values)
thstock   <- time * 0  # stock (kg), theoretical value
thoutflow <- time * 0  # outflow rate (kg/s), theoretical value

##### (4) computing by finite difference method
stock[1] <- stockinit
for(k in 1:length(time)){
  inflow[k]  <- inflowc
  outflow[k] <- stock[k] * outcoef
  if(k < length(time)) {
    stock[k+1] <- stock[k] + deltat * (inflow[k] - outflow[k])
  }
}

##### (5) computing theoretical values
thstockfinal <- inflowc / outcoef
for(k in 1:length(time)){
  thstock[k]   <- thstockfinal + (stockinit - thstockfinal)*exp( -outcoef*time[k] )
  thoutflow[k] <- thstock[k] * outcoef
}

##### (6)  output 
##### (6a) time-series graph of stock
# plot(time, stock, tck=1)
# lines(time, thstock, col='green')

##### (6b) time series graph of flow
 plot( time, inflow, tck=1)
 points(time, outflow, pch=2)
 lines(time, thoutflow, col='green')

2015-06-02, 更新(最新) 2015-06-09
増田耕一 (MASUDA Kooiti)

[システムに関する議論の目次へ] [この科目の教材目次へ]