線形微分方程式系の直接解法(ボックスモデル)

　今日は屋内での汚染物質濃度を求めることを仮定した簡易ボックスモデルを例にとって定数係数の連立線形常微分方程式の数値解を求めるための１つの手法を紹介します。

　室内の濃度分布を計算するには簡易モデルではなくて熱対流などの気流を計算する２方程式モデルなどに移流拡散方程式を付加した非定常な数値流体方程式を解くという大げさな方法もあります。

　しかし,ここでは部屋の平均濃度だけを計算する,という簡易ボックスモデルを考察します。

　例えば容積がＶ(ｍ³)の部屋の中央に石油ストーブなどがあって,それからの一酸化炭素(CO)の排出強度がＱ(ｍ³/s)であるとします。

　もしも部屋が密閉されていれば,時間ｔが経った後に,その一酸化炭素の平均濃度をＣとすると,Ｃ＝Ｑｔ/Ｖになると考えていいでしょう。

　しかし一般に部屋には隙間があって換気されるわけで,通常１時間に部屋全体の空気が入れ替わる回数を換気率という量で表わします。

　例えば容積Ｖの部屋の換気率がｎなら,１時間の後には総量ｎＶの空気が排出され,物質量の保存則によって同時にｎＶの空気が外部から入ってくるわけです。

　それ故,１秒当たりの空気の流出入量をｋ(ｍ³/s)とするとｋ＝ｎＶ/3600 ですね。

そこで,各時刻の部屋の中の一酸化炭素の平均濃度をＣ(ｔ)で表わすとこの部屋の中の単位体積当たりのその量は各時刻にＣ(ｔ)(ｍ³/ｍ³)なので時刻Δｔの間に外気に流出する一酸化炭素の量はｋＣ(ｔ)Δｔです。

一方,外気の一酸化炭素濃度(一定:constant)をＣ₀とすると,

流入量はｋＣ₀Δｔですから,部屋内のその量の時刻変化は,

Ｖ(Ｃ(ｔ＋Δｔ)－Ｃ(ｔ))＝ｋ[Ｃ₀－Ｃ(ｔ)]Δｔなる式で

与えられます。

　さらに排出強度をＱとすると,右辺には発生量ＱΔｔが加わりますから,

結局,一酸化炭素濃度に対する微分方程式は,

ｄＣ/ｄｔ＝ａ(Ｃ₀－Ｃ)＋Ｑ/Ｖ　と表わすことができます。

　ここで,ａ≡ｋ/Ｖです。

　Ｂ＝ａＣ₀＋Ｑ/ＶとおけばｄＣ/ｄｔ＝－ａＣ＋Ｂとなりますが,

　この微分方程式は簡単に解けて,

　Ｃ(ｔ)＝(Ｃ( 0 )－ａ^-1Ｂ)exp(－ａｔ)＋ａ^-1Ｂ　

　なる解が得られます。

こうしたモデルを室内物質濃度の簡易ボックスモデルと言います。

　これを一般化して各部屋の容積がＶ_i(ｉ＝1,2,..,ｍ)のｍ個の部屋を有する住宅があって,これらの部屋の平均濃度が(ｍ＋1)次元の列ベクトルでＣ(ｔ)＝^t(Ｃ₀,Ｃ₁,Ｃ₂,..,Ｃ_m)で表わされるとします。

　ここで,Ｃ_kはｋ番目の部屋における物質濃度です。

　特に,Ｃ₀は屋外の平均濃度値を示しています。

そして各部屋には排出強度Ｑ_iの排出源があり"部屋iと部屋ｊの間の換気量＝先の方程式のｋに相当する値"をｋ_ij＝ｋ_ji(i≠ｊ;ｉ,ｊ＝0,1,2,..ｍ)とします。

すると,このボックスモデルの微分方程式は,

ｄＣ_i/ｄｔ＝∑_j=0^mａ_ij(Ｃ_j－Ｃ_i)＋Ｑ_i/Ｖ_i(ｉ＝0,1,2,..,ｍ)

となります。ただしａ_ij≡(ｋ_ij/Ｖ_i)です。

　さらに,(排出強度/容積)の列ベクトルを,

　Ｂ≡^t(Ｑ₀/Ｖ₀,Ｑ₁/Ｖ₁Ｑ₂/Ｖ₂,..,Ｑ_m/Ｖ_m)で定義します。

　ただしＱ₀/Ｖ₀＝0 としておきます。

　また,λ_i≡∑_j=0^mａ_ijと定義し,行列Ａを成分;(ａ_ij－λ_iδ_ij)によって定義すれば,濃度を求めるボックスモデルの微分方程式:(ｍ＋1)変数の定数係数連立線形非同次微分方程式は,

　(ｍ＋1)次元ベクトル空間での線形非同次微分方程式として,

　ｄＣ/ｄｔ＝－ＡＣ＋Ｂ　と表わされます。 　

　この方程式の形式的な解は簡単に求めることができて,

　初期値Ｃ(0)に対し,Ｃ(ｔ)＝{Ｃ(0)－Ａ^-1Ｂ}exp(－Ａｔ)＋Ａ^-1Ｂ

　と書くことができます。

　さらに,濃度をこうした形式解ではなく実際に具体的に求めるには,

　従来の慣用的な方法であれば,Ａの最大(ｍ＋1)個の固有値:λ_i(ｉ＝0,1,2,．．ｍ)を全て求め,Ｃ(ｔ)の各成分を具体的にexp(－λ_iｔ)の

１次結合で表わします。

　しかし,高々10部屋程度の住宅で行列Ａが具体的にわかっているときには,固有値を求めるなどという面倒な手続きを省略して,コンピュータで無理矢理,力技で行列の指数関数を計算してしまう,という方法があります。

　形式解:Ｃ(ｔ)＝{Ｃ(0)－Ａ^-1Ｂ}exp(－Ａｔ)＋Ａ^-1Ｂで,

　逆行列:Ａ^-1はＡからコンピュータで掃き出し法などで計算可能ですし,

　exp(－Ａｔ)＝∑_ν=0^∞(－Ａｔ)^ν/ν!であり,コンピュータで行列のベキ乗(power):(－Ａｔ)^νを具体的に計算することも可能です。

　そこで,左辺のexp(－Ａｔ)を右辺の∑_ν=0^∞のν＝0,1,2,..の最初の数項で近似評価することによって力技で数値的にＣ(ｔ)の近似解を求めることができるわけです。

　逆に,この近似解と物質濃度の直接測定実験の結果から非線形最小二乗法などによって換気率パラメータ:ｋ_ijを推定することもできます。

http://fphys.nifty.com/(ニフティ「物理フォーラム」サブマネージャー）                                       TOSHI 

http://blog.with2.net/link.php?269343(ブログ・ランキングの投票）↑ここをクリックすると投票したことになります。