地球環境を考えるうえで、 物質・エネルギー循環系と フィードバック・システムの特徴を兼ね備えたシステムを考えることが重要である。 その簡単な例として、ここでは、「ふろおけモデル」を考える。
【ふろおけは、Meadows (2008)の入門書で、最初の例として出てくるが、 そこでは流出量が一定になっており、フィードバック・システムの例ではない (フィードバックは別の例を使って説明されている。) 他方、戸川 (1976)の本の、 計算プログラムをともなう例のうち最初に出てくるもの(第3章3.1節)が、 「水時計」となっているが、ここでいうふろおけモデルにあたるものである。 その本文では、水平断面積(下のわたしの記述での A)が高さによって違いうる場合について述べているが、 最初のプログラム例では、水平断面積も一定としている。】
入れものがあって、水がはいっているとする。 具体的に想像するものは、たとえば台所の流しでもよいのだが、 ここでは「ふろおけ」としておく。 (ふろがまとつながっているといった、ふろおけ特有の構造を考えるわけではない。) その形は、水平断面の形が一定のものを考える。 (ふろおけといえば、下のほうがつぼまったものを想像する人もいるかもしれないが、 ここではそうでないことに注意してほしい。) 次に示す図は、鉛直断面と思ってほしい。 底に小さなあながあいているとし、このあなには、栓をしないとする。 (住宅のふろでは考えにくいが、銭湯などではありうる状況だろう。)
このシステムの状態は、長さの次元をもつ「水位」h というひとつの量で表現できる。 たまっている水の体積は、水位に、一定と仮定した水平断面積 Aをかけたものである。 水の密度ρも一定と仮定するので、たまっている水の質量 S = ρ A h も、水位に比例する。 (記号 S は「たまり」にあたる stock のつもりでつけた。)
この、たまっている水の質量 S が、質量保存の法則に従って変化する。 流入量Fin、流出量Foutは、単位時間あたりの質量の流れとして考える。 流入量は、系外の条件によって決められるとする。 流出量は、系の状態つまり たまりの量 だけによって決められるとし、 たまりの量が多いほど多い(たまりの量の単調増加関数 Fout = Fout(S)である)とする[注1][注2]。
dS/dt = ( Fin - Fout(S) )
(1) まず、流入量がシステムの外の条件で決められた、 時間とともに変化しない一定値 Fin,1をとる場合を考えよう。 流出量がこれに等しくなる たまりの量 S1がある (Fout(S1) = Fin,1)としよう。 もし、現在のたまりの量 SがS1よりも高ければ、 流出量Fout(S)は流入量よりも大きいので、たまりの量は減るほうに変化する。 たまりの量がS1よりも少なければ、たまりの量はふえるほうに変化する。 たまりの量がS1ならば、その状態は変わらずに続く。 つまり、この系には、 定常入力 Fin,1を仮定すると、定常状態 S1がある (図の左の部分参照)。
(2) システムの外の条件が変わって、流入量がFin,2となり、そのまま持続したとする。 たまりの量は、新しい定常入力に対応する (Fout(S2) = Fin,2となるような) 定常状態 S2に近づくように変化するだろう。 (もし、流出量がたまりの量に比例するならば、近づきかたは、 新しい定常状態との差が指数関数(指数の係数が負)で減衰していくようなものになるだろう。)
(3) 出口の条件が変わって、たまりの量から流出量を決める関数 Fout(S)の定量的性質が変わったとする。 たとえば、出口が詰まりぎみになった場合の例として、 たまりの量 S から流出量を決める関数 Fout(S)に 1よりも小さい係数 b がかかった場合を考えよう。 流入量は一定値Fin,1のまま変わらないとする。 関数の変化の前に、たまりの量が定常状態にあったとすれば、 出口が詰まりぎみになると、流出量が減るから、水位は上がって、 前よりもたまりの量が大きい、新しい定常状態に向かうだろう(図の中央と右の部分参照)。
「ふろおけモデル」のうちでも 簡単な場合として、流入量(単位時間あたりの質量)Finは一定値であり、 流出量(単位時間あたりの質量)Foutが ふろおけにたまっている水の質量(以下、便宜上「たまりの量」と呼ぶ)に 比例する場合を考える。
式では次のように書ける。係数 c (便宜上「流出の係数」と呼ぶ)は定数である。
dS/dt = Fin - Fout ... (1)
Fout = c S ... (2)
(1), (2)をあわせてFoutをひとまず消去し、 Sは時間tの関数、Finとcは定数であることに注意すると、 このモデルは次のようなSに関する常微分方程式で表現できる。
dS/dt = Fin - c S ... (3)
ここで、仮に dS/dt = 0 という条件を置いてみると、それを満たす解(定常状態 steady state)がある。 添え字「steady」をつけて示すことにする。 Fout,steady = c Ssteady = Finであることから、 Ssteady = Fin / c である。 これを使って、(3)を S - Ssteady = S' についての式に書きかえると
dS'/dt = -c S' ... (4)
となり、その解は、指数関数 exp(x) = ex を使って
S'(t) = S'(0) exp( -c t ) ... (5)
のように書くことができる。ただしS'(0)は時刻 t = 0 でのS'の値である。Sについての式にもどすと
S(t) = Ssteady + (S(0) - Ssteady) exp( -c t ) ... (6)
となる。cは正の定数なので、(6)の右辺第2項は、時間 t が大きくなると 0 に近づく(減衰する)。 つまり、初期値S(0)がどんな値であっても、Sは時間とともにSsteadyに近づく。
このように、この例では、微分方程式が数式の形で解ける(「解析解が得られる」と表現されることが多い)。
モデルが少し複雑になると、解析解が得られなくなる。 その場合に有効なのは、数値計算によって近似解を求めることである。 有限の時間区間(時間差Δt)の間のSの変化 ΔS = S(t+Δt) - S(t) を考える。 dS/dtは ΔS/Δt の Δt→0 での極限だが、これを、極限をとらないΔS/Δtで近似する。 たとえば (1) の近似は次のようになる。
ΔS/Δt = Fin - Fout ... (7)
これを変形して、S(t)がわかっているときにS(t+Δt)を求める式が得られる。
S(t+Δt) = S(t) + ( Fin - Fout ) Δt ... (8)
いまFinは定数としたが、Foutは時間とともに変化する。 (8)式を使おうとすると、時間区間のうちのどの時刻でのFoutの値を使うかが問題になる。 ここでは、区間の初めの時刻 t での値で代表させることにする。
【微分を有限区間の差分で近似する方法にはいろいろなものがある。 気象学で応用される方法の概略は、別の教材ページ 「数値計算の教材: 時間発展問題の差分計算」で紹介した。 区間の中央で代表させる方法もあり、うまく使えればそのほうが精度がよいのだが、 数値があばれること(「計算不安定」という)をふせぐのにくふうがいるので、 ここでは、区間の初めで代表させる(専門用語で「前方差分」という)方法に限って紹介する。】
Δtが大きすぎると、差分による計算結果は微分方程式の解の近似として不適切になる。 Δtを小さくすると、計算量がふえて、計算機に対する負荷が重くなるほか、誤差の累積が問題になる場合もある。 適切なΔtを選ぶには、試行錯誤を必要とすることもあるが、 まず重要なのは、起こる現象の時間規模よりもやや細かくとることである。 たとえば、(6)式のような指数関数型の変化が起こるならば、Δtは 1/c よりも小さく(たとえばその5分の1くらいに)する。
この簡単な例を、R言語でプログラムにしてみる。 Rは統計的データ解析のための言語であり、 ここで行なうのは数値シミュレーションの基礎となる理論的計算なので、 Rの典型的な使いかたではないのだが、なんとか使うことができる。 Rの処理系は、2019年度から、板倉キャンパスのパソコン教室のWindowsマシンに はいっているし、 オープンソースソフトウェアとして提供され、 MS Windows, Mac, Linuxをはじめとする多くのOSの計算機上で動くので、 Rに慣れることは多くの進路の人にとって損にならないと思う。
なお、Rを含むほとんどのプログラム言語は、数値を扱うことはできるが、 単位のついた物理量を扱う機能をもっていない。 物理量の次元と単位についてはプログラミングの前の段階で人が注意して扱わなければならない。
教材用のプログラムを作るにあたって、次のような方針をとることにする。
【コメントやグラフに表示される語句に、日本語が含まれている。 (もし英語のほうがよい人がいれば英語版もつくるので、講師にたずねてほしい。) パソコン教室では OS が Windows 7 なので、文字コードは「シフトJIS」のほうを使ってほしい。 新しい OS では、文字コードが「UTF-8」のほうを使う必要があるかもしれない。】
R言語では各行の「#」という字からあとはコメントであり計算処理上は無視される。
R言語で「<-」は代入をあらわす記号である。右側の式の計算結果を左側の変数にいれる。
例文のプログラムを部分に分けて説明するために、コメントに(1)(2)のような番号を入れた。(式の番号ではない。)
流入量 Fin (変数名inflowc)と、流出の係数 c (変数名outcoef)を与える。 inflowc の値は 10 kg/分 とし、単位を kg/秒 に換算する。 outcoef は直観的に与えにくいので、 この例に限っては解析解がわかっていることを利用し、 定常状態となるたまりの量 Ssteady (変数名stockstd)の代表例を 50 kgと決め、 それが定常状態となる条件から計算された流出の係数の値(の代表例)を outcoefstd とする。 例文では代表例の値をそのまま outcoef に代入している。
差分の1ステップの時間間隔Δt (変数名deltat)を1分、 シミュレーションのステップ数(変数名 nstep)を20 (つまり継続時間を20分)とする。 (シミュレーションの考えかたのうえでは継続時間から考えて割り算でステップ数を決めるべきだが、 割り算の結果が整数にならないとやっかいなので、この手順にした。) また、初期(t = 0)でのストック量S(0) (変数名 stockinit) を 30 kgとしてみる。 (これをSsteadyと同じにするとつまらない結果になるので、 それと同じ桁で別の値を入れてみる。)
時間とともに変化する変数(および、この例では一定だが、別の問題では変化させたい、流入量inflow)に ついて、それぞれ、(nstep+1)個の記憶場所を確保する。 この手続きの書きかたは、プログラム言語によって違う。 なお、(3x)と注記した部分は、(5)の理論値(theoretical values)の計算のためのものなので、 ふつうのシミュレーションの場合は、なくてよい。
シミュレーションの本体。まずストック量の初期値を与えたあと、 for文によるループで、時間をΔtずつ進めながら、 nstep回、次の時間ステップでのストック量を計算する。
【この部分のプログラムは、R言語の「for(k in 1:length(time))」のループで、 kの値を1から順にひとつずつ length(time) [この例文では21] まで変えて それぞれ実行されることを前提としている。 もし並列に、あるいは順不同に実行される場合は 期待どおりの動作をしない。 その場合は確実に順に実行されるようにプログラムを変更する必要がある。】
この例題には、上に述べたように、解析解があるので、それも計算してみる。 (解析解があれば差分計算する必要はないのだが、両者の結果を比較することによって 差分計算の方法の精度を評価することができ、 解析解がない問題についての差分計算の結果を評価する際に役だつ。)
結果をグラフにする。ただし、画面表示ではなく、PNG画像ファイルで書き出すことにした。 流入量・流出量の時系列のグラフと、 たまりの量の時系列のグラフとを、2つの別々の画像ファイルにする。 差分による計算結果を記号(まるや三角)、理論値を折れ線で表示することにする。
流入量と、流出の係数の値を、Rのplot機能の「subtitle」として図の下に書くようにした。 現状の例文では、それと横軸のラベルとの間がせますぎる。 対策がわかったら例文を改訂するつもりだが、ひとまずこのまま使ってみてほしい。
数値を数字で表示する手続きも用意したほうがよいのだが、 プログラム実行直後に、Rの処理系内のコマンド行として (コンソール画面で、プロンプト「>」のあとに) 「stock」や「outflow」などの変数名を単独で打てば、 その変数の値が表示される (配列の場合には列挙される)ので、 この例題の結果を確認することに限ったとりあえずの目的にはそれで充分かと思う。
自分の実習用のディレクトリ(フォルダ)を用意したうえで、 そこに、プログラム例 huro_ex1.R をダウンロードして、保存してほしい。
(Windows 日本語版の場合) R処理系を起動し、 「ファイル」の「ディレクトリの変更」で自分の実習用のディレクトリを選択したあとで、 「ファイル」の「スクリプトを開く」で、Rプログラムを選択し、 「編集」の「全て実行」で実行するのがよいだろう。 プログラムの編集は、Rの中の編集画面でして、 結果を「ファイル」から「保存」または「別名で保存」するのがよいだろう。 (Rの外のテキストエディタ (Windowsの「メモ帳」を含む)で編集することもできる。 その編集結果を使うには、「スクリプトを開く」で開きなおす。)
そのまま実行するほかに、(1), (2)の部分を次のように変更して結果を比較し、特徴をのべてほしい。
計算される内容を変更したら、(6a) (6b) で設定している画像ファイル名の 「ex1」の部分をそれぞれ変えて、結果どうしを区別できるようにしてほしい。
プログラムのどこを変えたら結果がどう変わったかを確認するには、 プログラムの複数のバージョンを比較できる形で持つ必要がある。 エディタで単純に保存すると上書きになってしまうので、「別名で保存」を使う、 あるいは変更前のプログラムをコピーして名まえを変えておく、などのくふうを、各自してほしい。 (バージョン管理を支援するソフトウェアを使うとよいのかもしれないが、わたしは(まだ)使っていない。)
プログラム言語 Python (パイソン)で同じ計算をする例文を追加した。
Windowsでは、PythonのIDLEというエディタに読みこんで、 (操作メニューが英語ならば)「Run」から「Run module」で実行できるはず。
文字コードは UTF-8である。
配列変数の要素番号が Pythonでは 0 からはじまる (R では 1 からはじまる)。 また、配列変数に要素を追加するには append 手続きをつかっている。 そこで追加された要素の番号を、「#」のあとのコメントでしめした。
残念ながら、結果を出力するところがまだしっかりできていない。 Pythonにはグラフ表示のライブラリプログラムがあるのだが、 板倉のパソコン教室には導入されておらず、計算機システム入れかえまでは導入はむずかしいそうだ。 そこで、数値を数字で書き出すことにしたのだが、 わたしがまだ Python でファイルを作成してそこに文字列を書きこむ方法をよく習得していないので、 配列変数のなかみをそのまま「シェル」画面に書きだしてしまっている。 結果を使うには、シェル画面から、テキストエディタなり表計算のワークシートなりへ、 コピー・ペーストする必要がある。 今後、ファイルに書式をととのえて書き出すように改良したいが、 いつまでにすると約束するのはむずかしい。
今回の問題の計算は、表計算ソフトウェア(Microsoft ExcelやLibreOffice Calc)のワークシートでも可能である。 [別の授業の教材「ふろおけモデルの簡単な例、表計算ソフトウェアによる計算シート例」]で、 ほぼ同じ問題をあつかっている (数値は同じでない)。 しかし、ワークシートのセルに計算式を入れる形をとるので、 計算式がはいっているのか数値がはいっているのかが見わけにくく、 さらに、計算式をセル間でコピーする際に、「相対参照」(式の位置がずれたのにあわせて参照さきもずれる)と 「絶対参照」(式の位置がずれても参照さきのセルはかわらない)とを使いわける必要があるので、 まちがえやすい。 (実際、昨年度わたしが作った別の問題のワークシートにまちがいが残っていて、受講生にめいわくをかけてしまった。) 時間発展型のシミュレーションを表計算ソフトウェアでやることは、おすすめしない。