[2018-06-25 暫定]
「ふろおけモデル」とその例題については、次のように別ページに書きました。
例題の説明は、Rか表計算ソフトウェアかどちらかのページを見てください。
ここでは、 (東洋大学板倉キャンパスのパソコン教室にあるWindowsマシンで動かすことを想定して) Python言語によるプログラム作成を試みました。 まだ暫定的なものです。
プログラムの演算内容は R 言語版と基本的に同じですが、 グラフ表示をやめ、配列変数を使うのをやめ、 ループ内で、差分計算の結果と、理論値(指数関数による計算)の両方を計算するようにしました。
実行方法は、プログラムをIDLEのFile - Openから読みこみ、 必要に応じて編集してから、 IDLEのRun - Run Module で実行します。 結果は Python Shell の画面に出るので、 別にエディタ(「メモ帳」など)を開いて、cut and paste します。
出力は、左の列から順に、 1:時間ステップ番号、2:時刻、3:たまっている水の質量(差分計算)、4:流入量、 5:流出量(差分計算)、6:たまっている水の質量(理論値)、7:流出量(理論値)。
(結果の cut and paste は不便なので、 ファイル出力に変えようと思っていますが、いつまでに改訂できるか未定です。)
# bathtubex1.py
# A bath-tub model
# A simple case: outflow is proportional to stock
# Python program 2018-06-25 MASUDA Kooiti
##### (0) using libraries
from math import exp # used for theoretical values only
##### (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 some variables for theoretical values)
thstockfinal = inflowc / outcoef
##### (4) main computation
##### ##### initial condition
stock = stockinit
#####
for k in range(1, nstep+1):
##### ##### computing by finite difference method (part 1)
time = k*deltat
inflow = inflowc
outflow = stock * outcoef
##### ##### computing theoretical values
thstock = thstockfinal + (stockinit - thstockfinal)*exp( -outcoef*time )
thoutflow = thstock * outcoef
##### ##### output
print(k, time, stock, inflow, outflow, thstock, thoutflow)
##### ##### computing by finite difference method (part 2)
if(k < nstep):
stock = stock + deltat * (inflow - outflow)
# bathtubex1s.py
# A bath-tub model
# A simple case: outflow is proportional to stock
# Python program 2018-06-25 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)
##### (4) main computation
##### ##### initial condition
stock = stockinit
#####
for k in range(1, nstep+1):
##### ##### computing by finite difference method (part 1)
time = k*deltat
inflow = inflowc
outflow = stock * outcoef
##### ##### output
print(k, time, stock, inflow, outflow)
##### ##### computing by finite difference method (part 2)
if(k < nstep):
stock = stock + deltat * (inflow - outflow)