地球環境を考えるうえで、 物質・エネルギー循環系と フィードバック・システムの特徴を兼ね備えたシステムを考えることが重要である。 その簡単な例として、ここでは、「ふろおけモデル」を考える。
【ふろおけは、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くらいに)する。
この簡単な例を、Python (パイソン) 言語でプログラムにしてみる。 Python の処理系は、 オープンソースソフトウェアとして提供され、 Linux, MS Windows, MacOS をはじめとする多くのOSの計算機上で動く。
なお、Python を含むほとんどのプログラム言語は、数値を扱うことはできるが、 単位のついた物理量を扱う機能をもっていない。 物理量の次元と単位についてはプログラミングの前の段階で人が注意して扱わなければならない。
教材用のプログラムを作るにあたって、次のような方針をとることにする。
各行の「#」という字からあとはコメントであり計算処理上は無視される。
Python で「=」は代入をあらわす記号である。右側の式の計算結果を左側の変数にいれる。 なお、「==」は「等しい」という論理演算子である。
例文のプログラムを部分に分けて説明するために、コメントに(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)に ついて、それぞれ、個数がふえることのできる記憶場所を確保する。 この手続きの書きかたは、プログラム言語によって違う。 (Python では、空 (から) のリストをつくっておいて順に要素を追加することが 可能なので、ここではその方法をとることにした。) なお、(3x) と注記した部分は、(5) の理論値 (theoretical values) の計算のためのものなので、 ふつうのシミュレーションの場合は、なくてよい。
シミュレーションの本体。まずストック量の初期値を与えたあと、 for文によるループで、時間をΔtずつ進めながら、 nstep回、次の時間ステップでのストック量を計算する。
【この部分のプログラムは、Python 言語の 「for k in range(nstep+1)」のループで、 kの値を 0 から順にひとつずつ nstep [この例文では20] まで変えて それぞれ実行されることを前提としている。 もし並列に、あるいは順不同に実行される場合は 期待どおりの動作をしない。 その場合は確実に順に実行されるようにプログラムを変更する必要がある。】
この例題には、上に述べたように、解析解があるので、それも計算してみる。 (解析解があれば差分計算する必要はないのだが、両者の結果を比較することによって 差分計算の方法の精度を評価することができ、 解析解がない問題についての差分計算の結果を評価する際に役だつ。)
結果をグラフにする。 流入量・流出量の時系列のグラフと、 たまりの量の時系列のグラフとを、2つの別々の図にする。 差分による計算結果を記号 (まるや三角)、理論値を折れ線で表示する。 流入量と、流出の係数の値を、図の表題の2行めに書く。
作図の部分の、「#画面表示」と注記した「plt.show()」を生かすと、図は画面に出る。 「#画像ファイル」と注記した複数の行を生かすと、図は画像ファイルに書きこまれる。
数値を数字で表示する手続きも用意したほうがよいのだが、まだしていない。
Pandas のデータフレームというデータ構造をつかい、 書き出したいリストの内容をデータフレームにまとめて、 to_csv で CSVファイルに書き出すのがよさそうだと考えている。 【そのばあい、シミュレーションをおこなうためのリスト変数を、 データフレームの内わけと兼用するか、別に確保して内容をコピーする形にするかで、 まよっている。】
プログラムを対話型環境 (spyder や IDLE) で実行しているのならば、 プログラム実行直後に、その変数の値を表示させて見る方法はある。 (残念ながらわたしはその説明を書けるほど習熟していない。)
自分の実習用のディレクトリ (フォルダ) を用意したうえで、 そこに、プログラム例 huro1.py をダウンロードして、保存してほしい。
そのまま実行するほかに、(1), (2)の部分を次のように変更して結果を比較し、特徴をのべてほしい。
計算される内容を変更したら、(6a) (6b) で設定している画像ファイル名の 「huro1」の部分に文字列を追加するなどして、結果どうしを区別できるようにしてほしい。
プログラムのどこを変えたら結果がどう変わったかを確認するには、 プログラムの複数のバージョンを比較できる形で持つ必要がある。 エディタで単純に保存すると上書きになってしまうので、「別名で保存」を使う、 あるいは変更前のプログラムをコピーして名まえを変えておく、などのくふうを、各自してほしい。 (バージョン管理を支援するソフトウェアを使うとよいのかもしれないが、わたしは (まだ) 使っていない。)