準線形計画法(ネットワーク輸送問題)(３)

　さて,本題である私のオリジナル:準線型輸送問題のアルゴリズムに入ります。

３.流量値依存の準線型コスト関数を持つ輸送問題を解くアルゴリズム 

２では,各枝ｅ_mの単位コストｃ_mが流量ｘ_mの値に関係なく,一定であるようなごく一般的な線型輸送問題を扱った。

本節では,ｃ_mの値がｘ_mと共に単調増加し,コスト関数が下に凸の折れ線となるような輸送問題が解けるようにアルゴリズムを拡張する。 

すなわち,各枝は容量制限Ｌ_m≦ｘ_m≦Ｕ_m (Ｌ_m,Ｕ_mは±∞,正,負,0のいずれでも良い)を有するが,この他にｉ＝－ｊ_m,..,0,..,ｋ_mの各値について,ｔ_m,iが存在してｔ_m,i-1≦ｘ_m≦ｔ_m,iの各ｘ_mに対して,単位コストｃ_m,iを持つ。なお,ｔ_m,0＝0 である。

ｅ_mに関わる総コストは,ｘ_m＞ 0 のときは,ｔ_m,lm-1≦ｘ_m≦ｔ_m,lmならばΣ_k=1^lm-1ｃ_m,k(ｔ_m,k－ｔ_m,k-1)＋ｃ_m,lm (ｘ_m－ｔ_m,lm-1)(3-1a)であり,ｘ_m＜ 0 のときは,ｔ_m,-lm-1≦ｘ_m≦ｔ_m,-lmならばΣ_k=0^lmｃ_m,-k(ｔ_m,-k-1－ｔ_m,-k)＋ｃ_m,-lm (ｘ_m－ｔ_m,-lm-1)(3-1ｂ)である。

まず,条件(1-1)と全ての容量制限を満足する１つの初期可能木解が得られたとする。

このとき.この可能木解で現在の各枝ｅ_mの流量ｘ_mの値によって暫定的な上界ｔ_m,kと,下界ｔ_m,k-1,および単位コストｃ_m,kを定めることができるので,とりあえずこれらの値を固定しておく。 

基底枝ｅ_m＝(ｉ,ｊ)に対しては,ｔ_m,k-1≦ｘ_m≦ｔ_m,kにおける単位コストｃ_m,kについて双対基底として,Π_j＋ｃ_m,k＝Π_i(3-2)となるように、{Π_n}(ｎ＝1,2,..,Ｎ)の組を計算する。

　次に,ｅ_m＝(ｉ,ｊ)が非基底のとき,ｘ_m＝ｔ_m.k-1ならばＤ_m＝Π_i－Π_j－ｃ_m,k ,ｘ_m＝ｔ_m.kならばＤ_m＝Π_j－Π_i＋ｃ_m,kとおく。 

固定した暫定的な上界,下界,単位コストに対して,全てのＤ_mがＤ_m≦0 となるまで,通常の主単体法の更新を繰り返す。

全てのＤ_mがゼロ以下となったときには,2-2で記述した論旨によって,暫定的な容量制限内に流量が留まっている限りにおいて,コストは最小になる。

つまり,全てのＤ_mがＤ_m≦0 なら,コスト関数においてｘ_mの±Δ_mの変化と共に変化し得る量Σｃ_m,kmｘ_mが(2-2)式と同様に変形されて,Δ_m＞0 の変化に対して必ず非減少となるからである。

そこで次に,得られた暫定的最小コスト解において,非基底枝ｅ_mに対してｘ_mが下界ｔ_m.km-1に等しい場合には,これを上界と見なして,新しいＤ_mの値としてのＤ_m^*をＤ_m^*＝Π_j－Π_i＋ｃ_m,km-1で与える。

Ｄ_m≦0 より－Ｄ_m＝Π_j－Π_i＋ｃ_m,km≧0 であり,単位コストｃ_m,kはｋに関して単調増加であるからＤ_m^*≦－Ｄ_mであるが,もしＤ_m^*＞0 ならｅ_mは基底に入る枝の候補となる。 

同様にｘ_mが上界ｔ_m.kmに等しい場合には,これを下界と見なして,新しいＤ_mの値としてのＤ_m^*をＤ_m^*＝Π_i－Π_j－ｃ_m,km+1で与える。

Ｄ_m≦0 より－Ｄ_m＝Π_i－Π_j－ｃ_m,km≧0 であり,単位コストｃ_m,kはｋに関して単調増加であるから,やはりＤ_m^*≦－Ｄ_mであるが,もしＤ_m^*＞0 ならｅ_mは基底に入る枝の候補となる。 

もし,以上のような操作で得られたＤ_m^*が,なおゼロ以下ならば既に得られている可能木解が求める最小コスト解である。

なぜなら,(3-1)式を,ｍ＝1からＭまで加えた総コストはΣ_m{Σ_k≦km-1ｃ_m,k(ｔ_m,k－ｔ_m,k-1)＋ｃ_m,kmｔ_m,km-1}＋Σ_ｅm∈ＵＤ_mｔ_m,km－Σ_ｅm∈ＬＤ_mｔ_m,km-1＋Σ_nΠ_nｒn－Σ_{￢ｅm∈Ｔ}Ｄ_mΔ_m(3-3)と書けるから,Ｄ_m≦0 によりΔ_m＞0 に対しては既に最小であり,上述の操作による非基底枝の境界を越えるΔ_m＜0 の方向へのｘ_mの変化に対しては,コスト関数を(3-3)式の形に変形した場合,ｘ_mの変化Δ_m^*＝－Δ_m＞0 についても,なお更新されたＤ_mの値Ｄ_m^*が全てゼロ以下となるからである。

しかも,上に示したように,境界を越える区間変更の操作に対して,Ｄ_mあるいはＤ_m^*は必ず非増加であるから,ある操作で全てのＤ_m^*がゼロ以下なら,続く操作でも常にＤ_m^*がゼロ以下という性質が保持される。

こうした説明からは,得られた解が最小コスト解ではなく,極小コスト解であることしか証明されないように見えるが,単位コスト関数が単調でありコスト関数が滑らかな変動をするため,Δ_mの値が大きく境界を越えた場合でも既述の計算式が若干の修正を除いて保持されることを考慮すると,全てのＤ_m^*がゼロ以下という結果が得られた場合,求める最小コスト解が得られたとして良いことがわかる。

一方,Ｄ_m^*＞0 の枝ｅ_mが発見されれば,それを基底に入る枝として新しい暫定的上界,下界,単位コストを設定し,基底から出る枝を選んで通常の主単体更新の手続きを繰り返せば良い。

実際の計算機におけるアルゴリズムも2-3に示した手順に暫定的上界,下界,単位コストを用いるという僅かな変更を行なうだけで容易に拡張することができる。

Ａ.付録　初期可能木解の便利な求め方

　輸送問題において、主要な課題である更新の手続きに関する方法を考察してきたが､実際のアルゴリズムでは(1-1)式と容量制限を満足する初期可能木解(初期極大木)を求めるのも重要な課題である。

本節では、V.Chvatalの方法を参考にして,一般的で簡単かつ便利な初期可能木解の発見方法を紹介する。

まず,もとのＮ個の節点に番号Ｎ＋1のダミーの人為節点を付加し,さらにもとの節点1,2,..,Ｎと,この人為節点を結ぶＮ本のダミーの人為枝ｅ_M+1,ｅ_M+2,..,ｅ_M+Nを付加する。

流量は,全て非負の値にするため,新しい流量変数をｘ_m'≡ｘ_m－Ｌ_m (ｍ＝1,2,..,Ｍ)(A-1)で与え,上界をＵ_m－Ｌ_mとする。

(1-1)より,ｒ_n＝Σ_ｅm∈Ａnｘ_m－Σ_ｅm∈Ｂnｘ_mであるから,新しい供給量ｒ_n'をｒ_n'＝Σ_ｅm∈Ａnｘ_m'－Σ_ｅm∈Ｂnｘ_m'＝ｒ_n－Σ_ｅm∈ＡnＬ_m＋Σ_ｅm∈ＢnＬ_m (ｎ＝1,2,..,Ｎ)(A-2)と変換して,ｒ_N+1'＝0 とおく。

　人為枝ｅ_M+iはｒ_i'＞0 のとき(ｉ,Ｎ＋1)で,ｒ_i'≦0 のとき(Ｎ＋1,ｉ)とする。

また,全ての枝ｅ_mの単位コストをｐ_m(ｍ＝1,2,..,Ｍ,Ｍ＋1,..,Ｍ＋Ｎ)として,ｅ_mが実枝ならｐ_mに 0 を人為枝なら 1 を与えることにより１つの線型輸送問題ができる。

　この最小コスト輸送問題を,もとの問題(主問題と呼ぶ。)に対して補助問題と呼ぶ。

　補助問題をまとめると,次のように表わせる。: 

　条件:Σ_ｅm∈Ａnｘ_m'－Σ_ｅm∈Ｂnｘ_m'＝ｒ_n'(ｎ＝1,2,..,Ｎ,Ｎ＋1) (A-3); 0≦ｘ_m'≦ Ｕ_m－Ｌ_m (ｍ＝1,2,..,Ｍ)(A-4a);0≦ｘ_m'(ｍ＝Ｍ＋1,..,Ｍ＋Ｎ)(A-4b)

　目的:ｚ＝Σ_m=1^M+Nｐ_mｘ_m'(A-5)を最小化する。

　ただし,ｐ_m＝0 (ｍ＝1,2,..,Ｍ)(A-6a)；ｐ_m＝1(ｍ＝Ｍ＋1,..,Ｍ＋Ｎ)(A-6b)である。

補助問題の初期可能木解はＮ個の基底枝を人為枝ｅ_M+1,ｅ_M+2,．．,ｅ_M+Nとして,それらの流量の初期値をｒ_i'＞0 ならｘ_M+i'＝ ｒ_i'(A-7a),ｒ_i'≦0 ならｘ_M+i'＝－ｒ_i'(A-7b)と設定し,Ｍ個の非基底枝ｅ_m(ｍ＝1,2,..,Ｍ)については,流量をｘ_m'＝0 (ｍ＝1,2,..,Ｍ)(A-8)と設定すれば得られる。

２で示した主単体法の更新を,この補助問題に反復適用して,得られた最小コストｚが正の数ならば,もとの主問題は不可能であって可能解は存在しないが,ｚがゼロとなる場合には主問題の可能解が存在する。

以下,ｚ＝0 の場合のみを考える。

補助問題の最小コスト解では,極大木に属する基底枝はＮ本であるから,少なくとも１本の基底枝は人為節点を端点に持つ人為枝である。

人為枝が１本だけの場合は,残りのＮ－1本の基底枝が元の主問題の基底枝となり,主問題の極大木を形成する可能木解がただちに得られる。

補助問題の基底枝の中に人為枝が２本以上ある場合には,ｚ＝0 よりそれらの人為枝の流量は全てゼロであるはずだから,人為枝の代わりに他のゼロ流量の非基底枝を加えてもｚ＝Σｐ_mｘ_m'に対する解としては,同値なものが得られる。

ここでの補助問題の極大木は,いくつかの実基底枝のみから成る部分木群を,人為節点Ｎ＋1を通る２本ずつの人為枝で結びつけたものとなるが,これらの部分木間の流量はゼロであるから,部分木同士を連結するために基底人為枝の代わりに,流量ゼロの適切な非基底枝をおきかえれば,実枝のみを基底とする主問題の極大木が得られるはずである。

ひとたび人為枝を含まない極大木Ｔが得られれば,補助問題の解の一部分である{ｘ_m'}(ｍ＝1,2,..,Ｍ)の各々に{Ｌ_m}の各々を加えることによって初期可能木解{ｘ_m}を求めることができる。

参考文献 

1.     A.I.Ali,R.V.Helgason and S.Lall:”Primal Simplex Network Codes:State-of-the-Art Implementation Technology.”Networks. Vol.8 (1978) pp315-339

2.     V.Chvatal(Ｖ.フバータル：阪田省二郎,藤野和雄,田口 東 訳)「線形計画法」(上)(下)(啓学出版)(1988)

3.     古林 隆:「線形計画法入門」(産業図書)(1980)

※     以上が1991年頃に仕事上で書いた大したこともない数値計算の１論文です。

実は,この頃私は友人と共同で理進ゼミという大学受験塾を始めて,この塾の講師兼経営者とある会社の嘱託社員の掛け持ちをしていました。

この塾の共同経営者で数学と化学担当講師は,本業が東京理科大事務官だったのですが,彼にこの論文とか「次元控除法によるポアソン方程式の解法」とか,その頃書いた論文を自慢げに見せたところ,「どこかのお金持ちたちが,最近もっとつまらない論文で理科大で博士号を取っている。これくらいの内容なら簡単に理学博士が取れるから推薦してあげよう。」と言われました。 

博士号なんて「足の裏の飯の粒」に過ぎないもので,これは「取らないと気持ち悪いが取っても大したことはない。」という例えですが,当時,私は修士号しか持っていないし,取っておけば将来これを使って詐欺師をやるわけではないが,何かの金儲けの種とかの役に立つかもしれない,という色気を出して彼に頼んでみました。

ただ申請が通っても38万円余りの手数料を取られるという話だったので本当かなという疑問も持ちました。

それで論文がこれだけでは少ないかなと思い,この２つに論文ではないけれど「風速と拡散係数が高さのベキの拡散方程式の解析解」という,かつて窒素酸化物総量規制マニュアルに載せる低煙源拡散式＝JEAモデル(環境庁モデル)を作った際に書いた報文と,それを中国人の研究員と役人に説明したときの英文の資料を加え,さらに1976年の修士論文「三重三元クォーク模型の束縛ポテンシャル」および当時指導教官であったＫ助教授と共に書いた「Quark Molecule」という昔の雑誌への投稿論文を加えて提出しました。

結局,教授会まで行ったけど,その最終審査で論文の数が少ない,という理由で却下されたらしいということだったので,ある意味でお金を取られないのでホッとしたという記憶があります。 

そもそも,問題は論文の数ではなくて中味だろうし,内容を読んでもいないし理解もしてないな,などと思いましたが,論文の内容にもそれほど自信がなかったので別にいいかなと思い直しました。

しかし,彼の話だとお金持ちならもっとつまらない論文でもかまわないらしい,ということだったし,私の内容なら間単に取れるということだった(随分話が違う)ので「私学というか,この世界も所詮はお金か」と改めて感じたものです。

もっとも,自慢じゃないけど私は免許とか資格とか"お上から頂くもの"は全て大嫌いなので,それを取ろうとしたこともほとんどなくて,そういうものは卒業証書以外にはありませんから,基本的には生活費の足しになるという理由以外では博士も不要です。

まあ,私も就職してからは,論文やレポートは専門とは程遠い,流体力学関係か数値計算関係のものばかり書いていました。

そういえば1976年の修士論文「三重三元クォーク模型の束縛ポテンシャル」や「Quark Molecule」も,心ならずも私の本当の専門である量子電磁力学(QED)や場の量子論(QFT)におけるカイラル・カレントに関わるWard-Takahashi identityから生ずるAdler-Jackew anomaly,すなわち,トライアングル・アノマリーとは全く異なる題材のものです。

これらは,当時,本当の専門では修士論文の題材が見つからず,仕方なく丁度その頃発見されたJ/ψ粒子(これは後にチャームの発見と同定されました。)がカラー８重項のグルオンの１つではないか,などと推測をした結果であり,もともとは量子力学のシッフの古いアイデアを流用したものです。

つまり,これらはユニタリ群のカシミア演算子の既約表現に対する固有値の関数としてカラー・クォークや反クォークによるＮ体共鳴としての素粒子の質量公式を導くことにより,現在の技術ではエキゾチック粒子がなぜ発見されないか？という理由を群論的に考察したもので,私にとっては専門違いなので不本意なものでした。

(エキゾチック粒子とは,クォークそれ自身やクォーク４体以上のハドロン,あるいはカラー１重項以外の粒子です。)　

なお,「Quark Molecule」については1976年6月の「素粒子論研究」に掲載されています。

まあ,これを読むことが可能な人には実名がばれてしまうけど別に秘密じゃないから,ばれてもいいです。

私が自分の本当の細かい専門であると自認していた分野は,典型的な例としてはπ⁰ → 2γの崩壊確率として観測される量である量子電磁力学の摂動の１次の発散をくりこんだ結果として現われるアノマリーです。

(もしもカレントがカイラルでなければファりィの定理(Furry)によって,１次のような奇数次の摂動の計算結果は常にゼロになります。)

もっとも,これがカラー・クォークに関する話題と全く無関係というわけではありません。

アノマリーの存在はカラー・クォークのカラーが３色しかなく,したがってカラー対称性がユニタリ群としてはＳＵ(3)に限られることの証拠になっています。

また,トフフト(t'Hooft)によって素粒子の理論がこうしたアノマリーからフリーになるためには,クォークのフレーバーとレプトンの世代数が完全に一致しなければならないことも示されています。

こうした重要な事実も,このアノマリーの存在から言えることです。

(2006年５/11の記事「波動関数の位相と電磁場」参照。。)

ところで,素粒子論とか宇宙論とか純理論的なことをやってきても,そうした知識や経験は,それ専門の研究所か大学ではなくて普通の科学計算のシミュレーションやアセスメントをやる程度の民間会社では全然役に立たないのじゃないか,と思われるかもしれませんが,実際にはしっかりと役に立ちました。

まあ,いろいろとあるのですが,典型例としては４次元時空の非ユークリッド幾何学が中心の一般相対性理論でさえ,ちゃんと利用できるということがあります。

つまり,山や谷などの起伏のあるでこぼこの地表面が境界であるような空間領域での気流などを解析する際に,２次元曲面を表わすためにメトリック・テンソルを使うのが非常に有効な手段となり,気体の運動を測地線の方程式のようなもので表現できることがわかります。

無駄話が長くなりましたが,この項目については終わりにします。

http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャーです。

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

http://homepage2.nifty.com/toshis-kaiga-auction/健康商品の店 「TRS健康ランド」

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

準線形計画法(ネットワーク輸送問題)(２)

続きです。線形輸送問題のアルゴリズムの１例,を具体的に述べます。これはプログラムそのものを記述したものなので,少々冗長で退屈なものです。

実際にプログラムを組む必要がないなら,必ずしも全てを追跡する必要はないと思います。まあこんな感じのものである,ということで眺めてみてください。

2-3.計算機における実際のアルゴリズム

　計算機で実際に輸送問題を解くアルゴリズムの例としてAli,Helgason,Kennington,Lallによるものを紹介しておく。

　まず,節点の１つ：１を根と呼び固定しておく。

節点ｉの深さｄ_iとは,iから単純経路によって根に到るまでの枝の個数である。ｄ₁＝0 と定義する。

節点ｉの先行点ｐ_iとは,ｉから根までの単純経路上にあって,深さが(ｄ_i－1)の節点のことである。ｐ₁＝0 と定義する。

節点ｊは,ｊから根までの間に節点ｉがあるときｉの継続点と呼ばれる。Ｕ_iをｉの全ての継続点の集合とする。

一方,ｆを節点集合Ｖ＝{1,2,..,Ｎ}からＶの上への１対１の写像であるとする。Ｕ_i≠φでｆ(ｊ)＝ｉのとき,{ｆ(ｋ):ｊ＋1≦ｋ≦ｊ＋|Ｕ_i|}＝Ｕ_i (ただし,|Ａ|は有限集合Ａの元の個数)が成立すれば,この写像ｆを１つの継続点配列と呼ぶ。

１つの継続点配列が与えられたとき,節点ｉの糸ｓ_iとは,この配列でｉに続く節点,つまりｆ(ｊ)＝ｉならｆ(ｊ＋1)のことである。配列の最後の糸は根１であると定義する。

ｉに根付く極大木の最後の節点をｎ_iとすると,ｎ_iは次の条件のうちの１つを満足しなければならない。：(ａ)ｓ_ni＝１(ｂ)ｐ_ｓniは節点ｉから根までの道の上にある。

※(注)つまり,ｓ_ni≠１ならｐ_sni≠0 であり,ｓ_niから根までの道の上でのｓ_niの先行点がｐ_ｓniなので,ｐ_ｓniはもちろんｓ_niから根までを結ぶ道の上にある。

　そしてｓ_niはｎ_iから根１までをつなぐ|Ｕ_ni|個の点のうちのｎ_i以外のどれか１つである。

　しかもｎ_iはｉに根付く極大木の最後の節点だからｓ_niはi→ｎ_i→１という道の上にあり,ｐ_ｓniはｓ_ni→１なる道の上にあるから,結局ｐ_ｓniはｉと根１を結ぶ道の上にあることになる。

次にｅ_klがｉとｐ_iを結ぶ枝であるとする。

ｅ_kl＝(ｐ_i,ｉ)のとき,ｍ_i＝ｋ_i ,ｅ_kl＝(ｉ,ｐ_i)のとき,ｍ_i＝－ｋ_ｉ と定義する。ｅ_kl上の流量はα_iで表わすことにする。

基底に入る枝をｅ_φ＝(ｕ,ｖ)とするとき,更新手順は以下の「アルゴリズム１」および「アルゴリズム２」に示す通りである。

(１)アルゴリズム１

① 初期化

ａ. γ₁とγ₂を初期化する。

ｘ_φ＝Ｌ_φならγ₁＝－1,γ₂＝1；

ｘ_φ＝Ｕ_φならγ₁＝1,γ₂＝－1とおく。

ｂ．Δを臨界値としておく。

Δ＝Ｕ_φ－Ｌ_φ とおく。

ｃ.δ₁とδ₂を初期化する。

δ₁＝ｕ,δ₂＝ｖとおく。

② min(ｄ_u,ｄ_v)より大きい深さラベルを持つｕからｖへの部分経路上の最大許容流量変動Δを決定する。

ａ.最大深さラベルを持つ節点から始める。

ｄ_u－ｄ_v＜0 ならｄ＝ｄ_v－ｄ_u ,μ＝2とおく。ｄ_u－ｄ_v＝0 なら③へ進む。ｄ_u－ｄ_v＞0 ならｄ＝ｄ_u－ｄ_v ,μ＝1とおく。

ｂ．δ,γと計数を設定する。

δ＝δ_μ,γ＝γ_μ,ｋ＝1とする。

ｃ.増加流か減少流か？

ｍ_δγ＞0 ならｅへ進む。

ｄ.増加流

Ｕ_|ｍδ|－α_δ＞Δなら,δ＝ｐ_δとしてｆに進む。さもなければ,Δ＝Ｕ_|ｍδ|－α_δ,λ＝δ,δ＝ｐ_δとしてｆに進む。

ｅ.減少流

α_δ－Ｌ_|ｍδ|＞Δなら,δ＝ｐ_δとする。さもなければ,Δ＝α_δ－Ｌ_|ｍδ|,λ＝δ,δ＝ｐ_δとする。

ｆ.δの深さ＝min(ｄ_u,ｄ_v)

ｋ＜ｄならｋ＝ｋ＋1とおきｃに進む。

ｇ. δ_μの再設定

δ_μ＝δとおく。

③ｕからｖへの残りの道の上での最大許容流量変動Δの決定

ａ.道の終わりか？

δ₁＝δ₂なら,ｗ＝δ₁として終了する。さもなければｋ＝1とおく。

ｂ.δのロード

δ＝δ_k とおく。

ｃ.増加流か減少流か？

ｍ_δγ_k＞0 ならｅへ進む。

ｄ.増加流

Ｕ_|ｍδ|－α_δ＞Δなら,δ_k＝ｐ_δとしてｆに進む。さもなければ,Δ＝Ｕ_|ｍδ|－α_δ,μ＝ｋ,λ＝δ,δ_k＝ｐ_δとしてｆに進む。

ｅ.減少流

α_δ－Ｌ_|ｍδ|＞Δなら,δ_k＝ｐ_δとする。さもなければ,Δ＝α_δ－Ｌ_|ｍδ|,μ＝ｋ,λ＝δ,δ_k＝ｐ_δとする。

ｆ.両側を調べたか？

ｋ＝1ならｋ＝2として,ｂへ進む。さもなければａへ進む。

「アルゴリズム１」の終了時には,最大許容流量変動Δ,出る枝ｅ_|ｍλ|が定まり,ｕからｖへの道はｕからｗの道とｖからｗの道によって与えられる。

これに続く更新のアルゴリズムは次の通りである。

(２)アルゴリズム２

① 全更新が必要か否か？

Δ≠Ｕ_φ－Ｌ_φなら③へ進む。

②更新流の再設定

ａ.閉路に対する初期設定

δ₁＝ｕ,δ₂＝ｖ,ｋ＝1とする。

ｂ.枝ｅ_φの流量更新

ｘ_φ＝Ｕ_φならｘ_φ＝Ｌ_φとする。さもなければｘ_φ＝Ｕ_φとする。

ｃ.パラメータの初期化

δ＝δ_k,γ＝γ_kΔとおく。

ｄ.節点かｗか？

δ＝ｗならｈへ進む。

ｅ.増加流か減少流か？

ｍ_δ＞0 ならβ＝1,ｍ_δ＜0 ならβ＝－1とおく。

ｆ.更新流の再設定

α_δ＝α_δ－γβとおきかえる。

ｇ.次の先行点

δ＝ｐ_δとしｄに進む。

ｈ.ｋ＝1ならｋ＝2とおきｃに進む。さもなければ終了する。

③初期化

ａ.双対変数の変更

σ＝Ｄ_φとおく。

ｂ.入る枝の更新流量は？

ｘ_φ＝Ｕ_φなら,Δ'＝Ｕ_φ－Δ,σ＝－σとする。さもなければ,Δ'＝Ｌ_φ＋Δとする。

C.ｑ,ｑ'の初期化

μ＝1ならｑ＝ｕ,ｑ'＝ｖ,ε＝－φ,σ＝－σとする。さもなければｑ＝ｖ,ｑ'＝ｕ,ε＝φとする。

ｄ.ｉ,ｊ,γの初期化

ｉ＝ｑ,ｊ＝－ｐ,γ＝γ_kΔとおく。

ｅ.λの先行点のセーブ

λ'＝ｐ_λとしておく。

ｆ.ｑ'の深さラベルをセーブ

Ｆ＝ｄ_q'＋1としておく。

④ｉにおける双対変数、流量、枝ラベルの更新

ａ.双対変数の更新

Π_i＝Π_i＋σとおく。

ｂ.流量のセーブ

α'＝α_iとしておく。

ｃ.流量の更新

α_i＝Δ'とする。

ｄ.流れの方向のセーブ

ｍ_i＞0 ならβ＝1, ｍ_i＜0 ならβ＝－1とおく。

ｅ.枝のセーブ

ｍ＝|ｍ_i|と枝番号をセーブする。

ｆ.枝ラベルの更新

ｍ_i＝εとする。

ｇ. ｉの深さｄ_{i ,}およびＦとｄ_iの差をセーブする。

Ｆ'＝ｄ_i,Ｌ＝Ｆ－ｄ_iとしておく。

ｈ. ｉの深さの更新

ｄ_i＝Ｆとする。

⑤ｉの最後の継続点を決定し,双対変数と深さラベルを更新する。

ａ. ｉから始める。

ｚ＝ｉ,ｘ＝ｓ_iとおく。

ｂ. ｘは１つの継続点か否か？

ｄ_x≦Ｆ'なら⑥へ進む。

ｃ.双対変数とｘの深さの更新

Π_x＝Π_x＋σ,ｄ_x＝ｄ_x＋Ｌと更新する。

ｄ.続く糸

ｚ＝ｘ,ｘ＝ｓ_xとしてｂへ進む。

⑥節点ｉに対して,その"逆糸＝親"を発見する。

ａ. ｊから始める。

ｒ＝ｊとおく。

ｂ. ｒは逆糸か？

ｓ_r＝iならば⑦へ進む。

ｃ.続く糸

ｒ＝ｓ_rとしてｂに進む。

⑦先行点順列と糸の更新

ａ.最終更新か？

ｉ＝λならｉに進む。

ｂ. ｊに対する更新流量の決定

Δ'＝α'－γβとおきかえる。

ｃ. ｊの枝ラベルを決める。

ε＝－βｍとする。

ｄ. ｒとｚの糸の更新

ｓ_r＝ｘ,ｓ_z＝ｊとおく。

ｅ. ｉをセーブする。

ｒ＝ｉとおく。

ｆ. ｉの更新

ｉ＝ｊとおく。

ｇ. ｊの更新

ｊ＝ｐ_jとおく。

ｈ. ｉの先行点の更新

ｐ_i＝ｒ,Ｆ＝Ｆ＋1として④へと進む。

ｉ. ｒとｚの糸の更新

ｓ_r＝ｘ,ｓ_z＝ｓ_q'とする。

ｊ. ｑ'の糸の更新

ｓ_q'＝ｑとする。

ｋ. ｑの先行点の更新

ｐ_q＝ｑ'とする。

⑧ｑ'からλ'までの道の上の流量の更新

ａを次のように変えてｂを削除することを除いて②と同じである。

ａ.閉路に対する設定

μ＝1ならδ₁＝λ',δ₂＝ｑ',ｋ＝1とおく。さもなければ,δ₁＝ｑ',δ₂＝λ',ｋ＝1とおく。

以上である。

※     ここまでが,既存の「線型コスト関数に対する輸送問題のアルゴリズム」の１例です。今回はここまでにして,次回は「準線型コスト関数」に対するオリジナルな内容を述べる予定です。

参考文献

1.     A.I.Ali,R.V.Helgason and S.Lall:”Primal Simplex Network Codes:State-of-the-Art Implementation Technology.”Networks. Vol.8 (1978) pp315-339

2.     V.Chvatal(Ｖ.フバータル：阪田省二郎,藤野和雄,田口 東 訳)「線形計画法」(上)(下)(啓学出版)(1988)

3.     古林 隆:「線形計画法入門」(産業図書)(1980)

http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャーです。

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

http://homepage2.nifty.com/toshis-kaiga-auction/健康商品の店 「TRS健康ランド」

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

準線形計画法(ネットワーク輸送問題)(１)

　書類の整理をしていたところ,昔1991年頃に２番目の会社での社員時代に書いた未発表の数値計算に関する私の論文が見つかったので,それを紹介したいと思います。

これは,当時仕事上で,ある水道局の上水道について圧力損失が最小となるようなネットワークを構成する方法を開発してプログラム化する必要があって,その過程で書いたものです。 

こうしたことについて,当時私は全くの素人だったので,まず,線型計画法の専門書を数冊読みました。

その後,結局,線型問題ではなくて,より実際に近い近似的に２次の問題である準線形な"輸送問題＝最小コスト問題"を解くための当時オリジナル？と思われる"方法＝アルゴリズム"を自力で思いつき,それの論文を書いて某大手の中央研究所に収めました。

ネットワーク自体については,学生時代に専門の量子電磁力学(QED)や場の量子論で摂動計算を行なう必要があり,ファインマン・ダイアグラムを用いて具体的にＳ行列要素の摂動の２次や４次のグラフの計算をしたり,また,そうしたネットワークのグラフの解析性や構造を論じたりしたことはありました。

特にファインマン(Feynman)やダイソン(Dyson)らのくり込み理論においてグラフのツリー構造や"閉路＝ループ"の構造と摂動における運動量積分の発散の次数の関係を求めることにより,くり込み可能性を考察する上で様々なネットワーク・グラフの構造を調べた経験があります。

また,独学ですが主に線形素子から成る直流や交流の電気回路の理論においても,回路グラフなどを学習したことがあります。

それらはこの理論のネットワークとも共通した話なので,この理論には比較的違和感がなく入れました。

当時,本論文の方法を用いて実際にFortranでプログラムを組みましたが,これ自体の内容をプログラムするのは,そもそもこの論文がアルゴリズムの手順を書いたものですから,決して難しいものではありませんでした。

しかし,そのプロセスとして,いろいろな場合を想定して場合分けをする必要性があり,そのために順列,あるいは組合わせの全てを尽くす,という問題が生じました。

これは通常の数学の問題としては高校でも習う単純なものですから,最初は安易に考えていました。

しかし,実際にそれをプログラムにするためにアルゴリズム化しようとしたとき,非常に複雑な問題であることに気付き,結局,このプログラムを完成させるために,本質的内容には関係ないその部分でとても苦労した覚えがあります。 

もちろん,あらゆる場合についてネットワークのコストを全て計算して最小値を求めるというのでは,何の工夫もなく(非)線型計画法など全然必要ないわけですが,(非)線型計画法をうまく利用しても,なおかつ,かなりの場合分けが必要でした。

そして,実際,この順列,あるいは組合わせを尽くす,というプロセスのためには膨大な計算時間を必要とするので,当時利用できた今よりメモリーが少なく,ディスクのサイズも小さくて,CPUによる演算スピードも遅いコンピュータでは,計算時間がかかり過ぎる,というネックがありました。

プログラムの使用メモリ－のサイズを小さくしたり,計算時間を短縮したりするために,さらにいろいろと工夫した覚えがあります。

その後まもなく,私はその会社もやめたので,以後そのプログラムがどうなったかは知りません。 

論文のほうは,原理的には,これをプログラミングすることが可能であることを述べたものです。

私のやる仕事の特徴なのですが,ある意味で"理屈だけでくその役にも立たない"と言われるものの部類に属しています。

つまり,アルゴリズムであるにも関わらず,実用的であるというよりも単に理論的な意味しか持たないかも知れません。 

　論文の正確な表題は「準線形コスト関数を最小化する特殊なネットワーク輸送問題を解決するアルゴリズム」です。以下に論文をそのまま写します。

                        Abstract

Resently,because of the development in linear planning methods and the utility of digital computers with high efficiency, the techniques for solving the problems of minimizing network flow costs have been successfully developed by many investigators.

Most of these studys is due to the primal simplex algorithms.

In this paper,we modify the ordinary primal simplex method to apply the minimal cost network flow problems with more complex cost fuctions that increases with network flow amount on each arc.

1.     序文 

ネットワーク輸送問題,あるいは最小コスト問題は一般的な線型計画問題とは別に,通常の方法で線型方程式を解くことなく解ける問題であるという意味で,線型計画法の特殊な１分野として多くの研究者によって独立に研究されてきている。

もちろん,原理的には一般的な線型計画問題の直接的な応用として解くこともできるが,変数が多いときには輸送問題特有の方法を用いた方が、計算の簡略化,短時間化等でははるかに効率的であって,時間的,コスト的に非常に有利である。 

輸送型の問題としては,単純にＭ個の積み出し地からＮ個の目的地へ製品を輸送するのに必要な最小コストを求めるという,ヒッチコック型の問題が有名である。

しかし,本論文では,需要,供給を与えるＮ個の節点と,これらを任意に連結する独立のコストを持ったＭ個の枝から成る全く一般的な問題を扱うことにする。

各枝ｅ_m＝(ｉ,ｊ)は,節点ｉから節点ｊへの向きを持った流れを持つとする。(ｍ＝1,2,..,Ｍで,ｉ,ｊは1,2,..,Ｎの選ばれた異なる２つの番号である。)

これらの枝の集まりの全てが,実際に輸送可能な全てのルートを表わしている。一方,各節点ｎ(ｎ＝1,2,..,Ｎ)は付随する量ｒ_nを持つ。

ｒ_n＞0 ならｎはｒ_nに等しい供給を与える供給点である。ｒ_n＜0 ならｎは－ｒ_nに等しい需要を与える需要点である。またｒ_n＝0 ならｎは中間点である。つまり,ｒ_nは節点ｎにおける代数的供給量である。

任意の枝ｅ_m＝(ｉ,ｊ)を通過する流れがあるとき,単位流量当たりのコストがｃ_mであるとする。

このとき,最小コスト輸送問題とは,与えられた需要,供給と各枝の容量制限に合致しながら,輸送コストの合計を最小とするために,どの枝を用いて,それら各枝の流量をどのように決定するか？ということを定める問題である。

これらを数式化すると次のようになる。 

各節点ｎに対し,ｎから出て行く枝の集合をＡ_n,ｎに入ってくる枝の集合をＢ_nとするとき,各節点ｎ(ｎ＝1,2,..,Ｎ)に対して,Σ_ｅm∈Anｘ_m＝Σ_ｅm∈Bnｘ_m＋ｒ_n(1-1)が成立しなければならない。

ここに,ｘ_mは枝ｅ_m上の流量である。ｘ_mに適当な大きさの制限がある場合には,その制限式も付加する必要がある。

これらの条件の下で,最小コスト問題はΣ_m=1^Mｃ_mｘ_m(1-2)が最小の値を取ることを要求する。

ここでは,大前提として総需要と総供給の収支のバランスがとれているとする。つまり,Σ_n=1^Nｒ_n＝0 (1-3)と仮定する。

そうでない場合は,適切な需給を持ったダミー終点と各節点からこのダミー終点へのコストゼロの枝群を加えることで,需給バランスのとれているケースに帰着できるので,以下Σ_n=1^Nｒ_n＝0　とする。

　以下では,まず既存の主単体法によって,単純な線型コスト関数を持ち,各枝の流量に上界と下界のある場合の輸送問題を解くアルゴリズムの１例を紹介する。

　次に,本論文の主旨であるところの各枝の流量値の絶対値に依存して増加する単位コスト値を持ち,近似的に流量の２次関数型をした準線型コスト関数を持つ輸送問題へとそのアルゴリズムを拡張する。 

2-1.木と可能木解

①節点ｕとｖの間の道とは,節点ｗ₁,ｗ₂,..,ｗ_kと枝ｅ₁,ｅ₂,..,ｅ_k-1から成るネットワークで,ｕ＝ｗ₁,ｖ＝ｗ_k かつ枝ｅ_iの両端点がｗ_i,ｗ_i+1であるようなものである。

②閉路とは,道に始点と終点を結ぶ枝を１つ加えてできるネットワークである。 

③全ての節点の間に道があるとき,ネットワークは連結であるという。 

④木とは,連結でかつ閉路を含まないネットワークのことである。 

⑤ネットワーク{1,2,..,Ｎ}の"全ての節点を含む木":Ｔを問題としているネットワークの１つの極大木という。

　極大木はＮ－1個の枝を持ち,木に属さない枝の流量を勝手に決めることにより,条件(1-1)を満足する{ｘ_m}の組が決まる。

⑥可能木解とは,極大木Ｔに対応する可能解である。

　つまり、枝ｅ_mに対しＬ_mをその流量の下界,Ｕ_mを流量の上界とするとき,ｅ_mがＴに属さないならｘ_m＝Ｌ_m,またはｘ_m＝Ｕ_mであって,全てのｎについて,条件(1-1)を満足する{ｘ_m}の組(ｍ＝1,2,..,Ｍ)をＴに対する可能木解という。

⑦可能木解が得られたとき,全ての枝のうちで可能木解に対応する極大木に含まれるものはＮ－1本であり,これらの枝を基底という。

　極大木に含まれない枝を非基底という。

2-2.双対変数Π_nの決定と基底に入る枝の選択

　求めるアルゴリズムに到る最初のステップは,各節点ｎに対する双対変数 Π_nとして全ての基底枝ｅ_m＝(ｉ,ｊ)に対し,Π_j＋ｃ_m＝Π_i(2-1)となるような変数:{Π_n}(ｎ＝1,2,..,Ｎ)の組を計算することである。

　Ｎ個の変数Π_nに対し,(2-1)はＮ－1個の方程式系なので,解に定数だけの任意性があるため,例えばΠ_Ｎ＝0 と仮定して極大木の枝に沿って次々にΠ_nを決めることができる。

　次に,ｅ_m＝(ｉ,ｊ)が非基底枝の場合(すなわちｅ_mがＴに属さない場合:￢ｅ_m∈Ｔ)には,ｘ_m＝Ｌ_mのときＤ_m＝Π_i－Π_j－ｃ_mとおき,これらの枝の集合をＬ,ｘ_m＝Ｕ_mのときＤ_m＝Π_j－Π_i＋ｃ_mとおき,これらの枝の集合をＵとする。

　ただしＬ_m＝Ｕ_mのときには,どちらか一方を適当に選ぶ。

　このとき,全てのＤ_mがＤ_m≦0 を満足するなら,現在の可能木解が求める最小コスト解の１つである。

　なぜなら,Σ_ｅm∈Anｘ_m＝Σ_ｅm∈Bnｘ_m＋ｒ_n (1-1)の条件を満たしながら,非基底についてｘ_m＝Ｌ_m＋Δ_m,またはｘ_m＝Ｕ_m－Δ_m (Δ_m＞0)とｘ_mを変化させたとき,(それによって基底のｘ_mは自動的に変化する。)

Σ_ｅm∈T(Π_j－Π_i＋ｃ_m)ｘ_m＝0 だから,総コストは,Σ_m=1^Mｃ_mｘ_m＝Σ_m=1^M(Π_j－Π_i＋ｃ_m)ｘ_mm＋Σ_m=1^M(Π_i－Π_j)ｘ_m＝Σ_ｅm∈T(Π_j－Π_i＋ｃ_m)ｘ_m＋Σ_￢ｅm∈T(Π_j－Π_i＋ｃ_m)ｘ_m＋Σ_iΠ_iΣ_ｅm∈Aiｘ_m－Σ_jΠ_jΣ_ｅm∈Bjｘ_m＝Σ_￢ｅm∈T(Π_j－Π_i＋ｃ_m)ｘ_m＋Σ_nΠ_n(Σ_ｅm∈Aiｘ_m－Σ_ｅm∈Bjｘ_m)＝Σ_ｅm∈UＤ_mＵ_m－Σ_ｅm∈LＤ_mＬ_m＋Σ_nΠ_nｒ_n－Σ_￢ｅm∈TＤ_mΔ_m (2-2)

　そこで,全てのＤ_mがＤ_m≦0 のとき,任意の非基底ｅ_m上のΔ_m＞0 の変化に対して,Σｃ_mｘ_mは現在の状態より必ず非減少となるため,現在の可能木解が１つの最小コスト解となる。

　したがって,逆にＤ_m＞0 の枝ｅ_mがあれば,容量制限を越えることなくｘ_mを変化させることができる場合には,その変化によってコストは必ず減少する。

　つまり,ゼロより大きいＤ_mを持つ枝は基底に入る枝の候補である。

　ところで,Ｄ_m＞0 に対応する基底に入る枝を１つ決めて極大木に付加すると,その時点で１つの閉路ができてしまう。

　基底に入る枝ｅ_mの流量ｘ_mの変化を±Δとすると,Δ＞0 が大きいほどコストの減少は大きいわけであるが,ｘ_mの±Δだけの変化に対して閉路内の基底枝の全てがそれぞれ±Δずつ変化するときに限って,(1-1)式のΣ_ｅm∈Anｘ_m＝Σ_ｅm∈Bnｘ_m＋ｒ_nが満足される。

閉路に属する全ての基底についての容量制限Ｌ_m≦ｘ_m≦Ｕ_m(ｅ_m∈Ｔ)を満足しつつ,最大の値となるようにΔの値を定めることができる。

それにつれて,ｘ_m＝Ｌ_m,またはｘ_m＝Ｕ_mとなって,基底から出て非基底となる１つの枝が得られる。

こうして新しい極大木と可能木解が得られる。

入る枝に対してＤ_m＞0 ,出る枝に対してＤ_m＝0 であるから,このプロセスでコストは減少する。

そこで今度は,再び双対変数 Π_nを更新しなければならない。この手順は次の通りである。

Ｔに入る枝をｅ_φ,Ｔから出る枝をｅ_ψとすると,新しい極大木はＴ＋ｅ_φ－ｅ_ψと書ける。

Ｔ＋ｅ_φ－ｅ_ψから枝ｅ_φを除いたネットワーク(Ｔ＋ｅ_φ－ｅ_ψ)－ｅ_φ＝Ｔ－ｅ_ψは,枝ｅ_φの端点を１つずつ含む２つの互いに素な木から成る。

ｅ_φ＝(ｕ,ｖ)のとき,節点ｕを含む木をＴ_u,ｖを含む木をＴ_vとする。

ｅ_φの単位コストがｃ_φなら,σ＝Π_v－Π_u＋ｃ_φとおいて既知のΠ_nから更新されたΠ_nの値Π_n^*を次式によって求める。

Π_n^*＝Π_n＋σ(ｎ∈Ｔ_u),Π_n^*＝Π_n (ｎ∈Ｔ_v)(2-3)である。

これらが更新の手続きである。

以下,全てのＤ_mがゼロ以下になるまで上の手順を反復すればよい。

※今日はここまでにして続きは後日にします。

参考文献

1.A.I.Ali,R,V.Helgason and H.S.Lall:"Primal Simplex Network Codes:State-of-the-Art Implementation Technology." Networks.Vol.8(1978) pp315-339

2.V.Chvatal(Ｖ.フバータル；阪田省二郎,藤野和雄,田口 東 訳)「線形計画法」(上)(下)(啓学出版)(1988)

3.古林 隆：「線形計画法入門」(産業図書)(1980)

http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャーです。

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

http://homepage2.nifty.com/toshis-kaiga-auction/健康商品の店 「TRS健康ランド」

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

次元控除法によるポアソン方程式の直接数値解法(補遺)

Poisson方程式:△φ＝∇²φ＝ρを,微分方程式のままで解析的に解くのであれば,Laplace演算子(Laplacian):△＝∇²の逆演算子△^-1を求めて解として,φ＝△^-1ρを求めることに帰着します。

△^-1は微分演算子の逆演算子なので,これは一般に積分演算子です。

つまり,△_xＤ(ｘ－ｙ)＝δ(ｘ－ｙ)となるGreen関数;Ｄ(ｘ)を求める

ことが,△Ｄ＝1＝δ^なる積分演算子:Ｄ^＝△^-1を求めることに相当

します。

Ｄ(ｘ)はFourier変換を利用したり,Poisson方程式の斉次方程式であるLaplace方程式:△φ＝0 の解を調べることから得ることができます。

特に,球対称なGreen関数は,Ｄ(ｘ)＝－1/(4πｒ)(ただしｒ＝|ｘ|)で

与えられます。

これを,ｘ,ｙ,ｚで有限回偏微分したものは,一般にはGreen関数ではないですが,原点ｒ＝0 を除けば,－1/(4πｒ)と同じく調和関数,つまりLaplace方程式の解です。

それらは,物理学的には,双極子,四重極子,...etc.という風に,多重極子と呼ばれ,数学的には球関数,あるいは球面調和関数と呼ばれます。

そうして得られたＤ(ｘ)を用いて,φ(ｘ)＝∫Ｄ(ｘ－ｙ)ρ(ｙ)ｄｙとおけば,φは確かに△φ＝∇²φ＝ρの解になります。

つまり,△Ｄ＝１＝δ^ですから,Ｄ^はＤ^＝△^-1という積分演算子で

あり,φ＝△^-1ρ＝Ｄ^ρによって解φを求めるのが通常の解法です。

ここで,Ｄ^というのは関数Ｄ(ｘ)を"演算子＝積分演算子"とみなすことを意味しています。

　Poisson方程式:△φ＝∇²φ＝ρを差分方程式化すると,差分化によって左辺のLaplace演算子:△＝∇²は広義の行列:Ａとなり,微分方程式△φ＝∇²φ＝ρは,結局,連立１次方程式:Ａφ＝ρに帰着します。

　この連立１次方程式の解;φ＝Ａ^-1ρを,ガウズザイデル法のような

"緩和法＝繰り返し計算法"で解くのではなく,

　複雑で丸め誤差などが累積しやすいという数値計算に特有の欠点はあるのですが,直接法である消去法など,ＱＲ分解などを利用して色々と工夫してできるだけ精度良く,簡単に直接的に求めようと試みたのが,

　本題の「次元控除法によるポアソン方程式の直接数値解法」の目指す

ところです。 

そもそも,どこかに１つの大きな源があるようなPoisson方程式の解

である,太陽系の重力場やCoulomb場のようなものを求めるのが目的

であれば,それを求めるのに数値計算に頼る必要はないです。

しかし,室内気流に対する速度ポテンシャルに相当するような圧力場や,あるいは地球温暖化の源の１つとなる温排水などの影響による海流の変化のソースポテンシャルのような,中途段階で現われる量は,

通常はかなり複雑な条件下での２次元,３次元のPoisson方程式の解なので,それを解く場合は数値計算に頼らざるを得ません。 

しかし,本題の方法は対象とする領域の形が矩形や直方体に限られているので,建物の室内のような形状についての計算には適していても,一般の種々の形状の境界を持つ領域に対しては普通に考えれば,適用できないと思われるので,汎用性が少ないと感じられるかもしれません。 

しかし,領域が矩形や直方体ではなくても例えば２次元の円板や３次元球体のように,少なくとも単連結な領域であれば,それらは適当な変数変換によって互いに矩形や直方体に変換可能です。

あるいは,逆変換を行うことによって矩形や直方体を球面や球体に変換することもできます。 

実際,私のアルゴリズムによって作成したプログラムが,Poisson方程式の解を求めるものとなっているかを検証するための１手段として,

３次元メッシュの空間領域を十分大きく取り,直方体での計算結果を,直方体から球に座標変換することによって,解φが確かに球対称なCoulomb場:φ(ｘ)＝－ρ/(4πｒ)に一致することを確かめました。

というわけで,少なくとも単連結な領域における３次元以下のPoisson方程式であれば,私の示したこの方法は数値解を求めるのに有効なものであろうと考えられます。

つまり,トポロジー的に同相であるなら,一方で可能な操作は他方でも可能であるという例になっています。

(球面は境界がないので矩形と同相ではありませんが,球面を２つに分けると半球面は矩形と同相です。)

http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャーです。

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

http://homepage2.nifty.com/toshis-kaiga-auction/「健康商品の店 タクザイ」

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

次元控除法によるポアソン方程式の直接数値解法(2)

前記事の続きとして２次元Poisson方程式の問題から３次元Poisson方程式への拡張について記述した私自身の論文の続きを掲載します。

３次元の直方体領域での,３次元Poisson方程式は,

(∂²φ/∂ｘ²)＋(∂²φ/∂ｙ²)＋(∂²φ/∂ｚ²)＝ρ(ｘ,ｙ,ｚ)

(ａ≦ｘ≦ｂ,ｃ≦ｙ≦ｄ,ｅ≦ｚ≦ｆ)...(3-1)で与えられます。

Neumman的境界条件は,(∂φ/∂ｘ)|_x=a＝ｆ⁽¹⁾(ｙ,ｚ),

(∂φ/∂ｘ)|_x=b＝ｆ⁽²⁾(ｙ,ｚ),(∂φ/∂ｙ)|_y=c＝ｇ⁽¹⁾(ｘ,ｚ),

(∂φ/∂ｙ)|_y=d＝ｇ⁽²⁾(ｘ,ｚ),(∂φ/∂ｚ)|_z=e＝ｈ⁽¹⁾(ｘ,ｙ),

(∂φ/∂ｚ)|_z=f＝ｈ⁽²⁾(ｘ,ｙ)　...(3-2)　です。

　直方体を,各辺がｈ_x,ｈ_y,ｈ_zであるようなＮ_x×Ｎ_y×Ｎ_z個の微小セル

に分割します。

  ｘ₀＝ａ,ｘ_Nx＝ｂ;ｘ_i－ｘ_i-1＝(ｂ－ａ)/Ｎ_x＝ｈ_x(ｉ＝1,2,..,Ｎ_x);

  ｙ₀＝ｃ,ｙ_Ny＝ｄ；ｙ_j－ｙ_j-1＝(ｄ－ｃ)/Ｎ_y＝ｈ_y(ｊ＝1,2,..,Ｎ_y);

  ｚ₀＝ｅ,ｚ_Nz＝ｆ；ｚ_k－ｚ_k-1＝(ｆ－ｅ)/Ｎ_z＝ｈ_z(ｋ＝1,2,..,Ｎ_z)

　２次元の場合と同様に(3-1)式を差分化すると,

　　[(φ_i+1,j,k－2φ_i,j,k＋φ_i-1,j,k)/ｈ_x²]

　＋[(φ_i,j+1,k－2φ_i,j,k＋φ_i,j-1,k)/ｈ_y²]

　＋[(φ_i,j,k+1－2φ_i,j,k＋φ_i,j,k-1)/ｈ_z²]＝ρ_i,j,k

(ｉ＝1,2,...,Ｎ_x;ｊ＝1,2,..,Ｎ_y;ｋ＝1,2,..,Ｎ_z).

　..(3-3)となります。

　φ_0,j,k,φ_Nx+1,j,k,φ_i,0,k,φ_i,Ny+1,k,φ_i,j,0,φ_i,j,Nz+1は境界条件

　(3-2)式を用いて２次元の場合の(2-4)式と同様に消去することができ

　ます。

(3-3)式を境界条件を考慮して行列×列ベクトル＝列ベクトルの表現

に変換すると,

Φ_k+1－2Φ_k＋Φ_k-1＋{(ｈ_z/ｈ_y)²Ａ~＋(ｈ_z/ｈ_x)²Ｂ~}Φ_k＝Ρ'_k

(ｋ＝2,3,..,Ｎ_z－1)...(3-4a),　

Φ₂－Φ₁＋{(ｈ_z/ｈ_y)²Ａ~＋(ｈ_z/ｈ_x)²Ｂ~}Φ₁＝Ρ'₁(3-4b),

－Φ_Nz＋Φ_Nz-1＋{(ｈ_z/ｈ_y)²Ａ~＋(ｈ_z/ｈ_x)²Ｂ~}φ_Nz＝Ρ'_Nz　

(3-4c)となります。

Φ_k≡^t(φ_1,k,φ_2,k,..,φ_Ny,k),

φ_j,k≡^t(φ_1,j,k,φ_2,j,k,..,φ_Nx,j,k),

Ρ'_k≡^t(ρ'_1,k,ρ'_2,k,..,ρ'_Ny,k)　

と定義します。

ρ'_j,k≡^t(ρ_1,j,k＋ｆ⁽¹⁾_j,k/ｈ_x＋ｇ⁽¹⁾_1,kδ_j1/ｈ_y－ｇ⁽²⁾_1,kδ_jNy/ｈ_y

＋ｈ⁽¹⁾_1,jδ_k1/ｈ_z－ｈ⁽²⁾_1,jδ_kNz/ｈ_z,ρ_2,j,k

＋ｇ⁽¹⁾_2,kδ_j1/ｈ_y－ｇ⁽²⁾_2,kδ_jNy/ｈ_y

＋ｈ⁽¹⁾_2,jδ_k1/ｈ_z－ｈ⁽²⁾_2,jδ_kNz/ｈ_z,..,ρ_Nx,j,k－ｆ⁽²⁾_j,k/ｈ_x

＋ｇ⁽¹⁾_Nx,kδ_j1/ｈ_y－ｇ⁽²⁾_Nx,kδ_jNy/ｈ_y＋ｈ⁽¹⁾_Nx,jδ_k1/ｈ_z

－ｈ⁽²⁾_Nx,jδ_kNz/ｈ_z)

(ｊ＝1,2,..,Ｎ_y;ｋ＝1,2,..,Ｎ_z)...(3-5)です。

　(Ｎ_x×Ｎ_y)次の正方行列Ａ~,Ｂ~はＮ_x次の単位行列Ｅと,

　(2-7)で与えたＮ_x次の三重対角行列Ａを用いて,次のように

　表わすことができます。

　Ａ~は第１行目が(－Ｅ,Ｅ,0,..,0),第ｋ行目(ｋ＝2,3,..,Ｎ_y－1)

　が(0,..,Ｅ,－2Ｅ,Ｅ,0,..,0),第Ｎ_x行目が(0,..,0,Ｅ,－Ｅ)で

　与えられるＮ_x次正方行列を成分とするＮ_y×Ｎ_yの

　三重対角行列です。　...(3-6)

　また,Ｂ~は対角成分が全てＮ_x次正方行列Ａでそれ以外の成分は

　ゼロ　のＮ_x次正方行列を成分とするＮ_y×Ｎ_y

の対角行列です。　(3-7)

　　ただし,Ｎ_x次正方行列Ａは次の通りです。

　Ａ~を対角化する(Ｎ_x×Ｎ_y)次の直交行列をＳ~_Nyとすると,

　^tＳ~_NyＡ~Ｓ~_NyはＥをＮ_x次の単位行列として,対角成分が,

　σ₁Ｅ,σ₂Ｅ,..,σ_NyＥの対角行列になります。　...(3-8)

　ここに,

　σ_m＝－4sin²{(ｍ－1)π/(2Ｎ_y)} (ｍ＝1,2,..,Ｎ_y)...(3-9)

　 (Ｓ~_Ny)_il,jm＝(1/Ｎ_y)^1/2δ_ij (ｍ＝1；ｌ＝1,2,..,Ｎ_y),

　 (Ｓ~_Ny)_il,jm＝(2/Ｎ_y)^1/2cos{(2ｌ－1)(ｍ－1)π/(2Ｎ_y)}δ_ij

  (ｍ＝2,3...,,Ｎ_y；l＝1,2,..,Ｎ_y) (ｉ,ｊ＝1,2,..,Ｎ_x)　.(3-10)

です。

　また,前記事の(2-9)によって,Λ_Ｎx＝^tＴ_ＮxＡＴ_Ｎx (Λ_Ｎxは対角成分

　がλ₁,λ₂,..,λ_Ｎxの対角行列)　...(3-11)であり,　

　λm＝－4sin²{(ｍ－1)π/(2Ｎ_x)} (ｍ＝1,2,..,Ｎ_x)です。

　Ｔ~_Ｎyを対角成分が全てＴ_ＮxのＮ_y次対角行列 (3-12)と考えれば,

　^tＴ~_NyＢ~Ｔ~_Nyは対角成分が全てΛ_ＮxのＮ_y次対角行列

　となります。.(3-13)

　(Ｔ~_Ny)_il,jm＝(1/Ｎ_x)^1/2δ_lm (ｊ＝1；ｉ＝1,2,..,Ｎ_x),

　(Ｔ~_Ny)_il,jm＝(2/Ｎ_x)^1/2cos{(2ｉ－1)(ｊ－1)π/(2Ｎ_x)}δ_lm

　(ｊ＝2,3,..,Ｎ_x；ｉ＝1,2,..,Ｎ_x) (ｌ,ｍ＝1,2,..,Ｎ_y)

　..(3-14)です。

　以上から,Ａ~の直交変換:^t(Ｓ~_NyＴ~_Ny)Ａ~Ｓ~_NyＴ~_Nyは,対角成分が

　σ₁Ｅ,σ₂Ｅ,..,σ_NyＥの対角行列になります。..(3-15),

　一方,Ｂ~の直交変換:^t(Ｓ~_NyＴ~_Ny)Ｂ~Ｓ~_NyＴ~_Nyは,対角成分が全て

　Λ_Ｎxの対角行列になります。..(3-16)

それ故,Ｑ~≡Ｓ~_NyＴ~_Nyと置けば,Ｑ~による直交変換によって,

Ａ~,Ｂ~は同時に対角化できて,

^tＱ~{(ｈ_z/ｈ_y)²Ａ~＋(ｈ_z/ｈ_x)²Ｂ~}Ｑ~は,対角成分が,

　[(ｈ_z/ｈ_y)²σ_mＥ＋(ｈ_z/ｈ_x)²Λ_Ｎx]

の対角行列となります。(3-17)

　この得られた対角行列をΩとし,Ｐ_k≡^tＱ~Φ_k,Ｒ_k≡^tＱ~Ρ'_k

(ｋ＝1,2,..,Ｎ_z)とすれば,(3-4)のPoisson方程式は,

　Ｐ_k+1－2Ｐ_k＋Ｐ_k-1＋ΩＰ_k＝Ｒ_k(ｋ＝2,3,..,Ｎ_z－1) (3-18a),

　Ｐ₂－Ｐ₁＋ΩＰ₁＝Ｒ₁ (3-18b),

　－Ｐ_Nz＋Ｐ_Nz-1＋ΩＰ_Nz＝Ｒ_Nz (3-18c),

　すなわち,　　

　Ｐ_i,j,k+1＋{(ｈ_z/ｈ_y)²σ_j＋(ｈ_z/ｈ_x)²λ_i－2}Ｐ_i,j,k＋Ｐ_i,j,k-1

＝Ｒ_i,j,k(ｋ＝2,3,．,Ｎ_z－1),

　Ｐ_i,j,2＋{(ｈ_z/ｈ_y)²σ_j＋(ｈ_z/ｈ_x)²λ_i －1}Ｐ_i,j,1

＝Ｒ_i,j,1,

　 {(ｈ_z/ｈ_y)²σ_j＋(ｈ_z/ｈ_x)²λ_i －1}Ｐ_i,j,Nz＋Ｐ_i,j,Nz-1

＝Ｒ_i,j,Nz (ｉ＝1,2,．,Ｎ_ｘ;ｊ＝1,2,．,Ｎ_y)　

　...(3-19) となります。

　これはｉ,ｊを固定したとき,Ｎ_z次の三重対角行列の係数を持った

　１次元の連立１次方程式系であり,係数行列は,i＝ｊ＝1のとき非

　正則,その他の場合は正則ですから,２次元の場合と同様,容易に

　解けます。

　１次元の連立１次方程式を解いてＰ_kが得られれば,

　それに対して,Φ_k＝Ｑ~Ｐ_kとすれば,最終的な解Φ_k,

　またはφ_i,j,kが得られます。

　なお,Ｎ_x×Ｎ_y次の直交行列Ｑ~＝Ｓ~_NyＴ~_Nyの陽な形を書き下すと

　(Ｑ~)_ij,mn＝Σ_k=1^NxΣ_l=1^Ny(Ｓ~_Ny)_ij,kl(Ｔ~_Ny)_kl,mn

＝_l(Ｔ_Nx)_i,m(Ｔ_Ny)_j,n となります。

　ただし,m(Ｔ_Ｎ)_lm＝(1/Ｎ)^1/2 (ｍ＝1),

　(Ｔ_Ｎ)_lm＝(2/Ｎ)^1/2cos{(2ｌ－1)(ｍ－1)π/(2Ｎ)}

　(ｍ＝2,3,..,Ｎ)です。...(3-20)

　※付録(Appendix):特殊な三重対角行列Ａを対角化する直交行列

　第１行目が(－1,1,0,..,0),第ｋ行目(ｋ＝2,3,..,Ｎ－1)が

　(0,..,1,－2,1,0,..,0),第Ｎ行目が(0,..,0,1,－1)で与えられる

　Ｎ次の三重対角行列をＡとしすると,

　既に示したように,固有値問題Ａｘ＝λｘ(ｘ≠0)の解として,

　Ｎ個の異なる固有値λ＝λ_m＝－4sin²{(ｍ－1)π/(2Ｎ)}

　 (ｍ＝1,2,..,Ｎ)が得られます。

　そして,固有値λに属する固有ベクトルｘは,

　θ＝θ_m＝(ｍ－1)π/Ｎ (λ＝λ_m＝－4sin²(θ_m/2) )として,

　ｘ＝^t(ｘ₁,ｘ₂,．．,ｘ_N)の成分ｘ_kが,

　ｘ_k＝ａ[exp(iｋθ)＋exp{－i(ｋ－1)θ}}

　＝2ａexp(－iθ/2)cos{(2ｋ－1)θ/2}なる形で与えられること

　も,既に前記事の(注２)で見た通りです。

　そこで,λ_mに属する固有ベクトルを改めてｘ_m＝^t(ｘ_1m,ｘ_2m,..,ｘ_Nm)

　とおけば,ｘ_lm＝ｃ_mcos{(2ｌ－1)(ｍ－1)π/(2Ｎ)}です。

　そして,固有値λ₁,λ₂,..,λ_N はすべて異なるので,

　ｍ≠ｎなら実内積(ｘ_mｘ_n)は直交します,つまり(ｘ_mｘ_n)＝^tｘ_mｘ_n＝0

　を満たします。

　一方,(ｘ_mｘ_m)＝^tｘ_mｘ_m＝ｃ_m²Σ_l=1^Ncos²{(2ｌ－1)(ｍ－1)π/(2Ｎ)}

　＝Ｎｃ_m²/2　(ｍ＝2,3,..,Ｎのとき),Ｎｃ₁² (ｍ＝1のとき)ですから,

　ｘ_lm＝(2/Ｎ)^1/2cos{(2ｌ－1)(ｍ－1)π/(2Ｎ)}(ｍ＝2,3,..,Ｎ),

　(1/Ｎ)^1/2 (ｍ＝1)とおけば,(ｘ_ｍｘ_n)＝^tｘ_mｘ_n＝δ_mnと正規直交化

　されます。

　そこで,Ｎ×Ｎ行列Ｔ_Nで成分が(Ｔ_N)_lm＝ｘ_lmとなるものを与えると,

　このＴ_N は１つの直交行列で,Ａは直交変換^tＴ_NＡＴ_N＝Λ(Λは対角

　成分がλ₁,λ₂,..,λ_Ｎの対角行列)によって対角化されます。

　(以上)

参考文献：1)Masako Ogura：”A direct solution of Poisson's equation by dimention reduction method.”Journal of the Meteorogical Soc. Japan. Vol.47. No.4(1969).pp319-322(日本気象学会誌、第47巻4号,1969) 2)矢嶋信夫、野木達男 著「発展方程式の数値解法」(岩波書店) (1977)他

http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャ　ー

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

http://homepage2.nifty.com/toshis-kaiga-auction/「健康商品の店 タクザイ」

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

次元控除法によるポアソン方程式の直接数値解法(1)

　引越し後の整理が完全には終わってなかったので,書類の整理をしていたところ,昔1987年前後にある会社の正社員時代に書いた未発表の数値計算に関する私の論文の草稿が見つかったので,それを紹介したいと思います。

当時,仕事上でストーブなどの暖房装置が原因で生じる室内汚染について,一酸化炭素や窒素酸化物などの濃度測定実験の結果に対して理論的な裏付けを与える必要がありました。

そのため,気流,温度,濃度の３次元の空間的な分布が時間が経つにつれて非定常的にどのように変化していくか？を見るという目的で具体的に流体の運動方程式を数値的に解くということになりました。

そこで,乱流部分に関しては２-方程式モデルを仮定して,時間を１次元差分に分割し,空間を３次元の格子(メッシュ)に区切ることによって非定常な差分方程式を設定して,それを数値的に解いたわけです。

その際,解法の過程の１部分,つまり中間段階として圧力に相当する量をPoisson(ポアソン)方程式の解として得る必要性が生じました。

そして,この論文はそうしたPoisson方程式を解く際の数値計算を効率的,かつ高い精度で行なうことを目的として,当時新たに考案した手法を紹介したものです。

　論文の正確な表題は,「境界条件がノイマン的である場合の次元控除法による３次元ポアソン方程式の直接数値解法」

(A Direct Numerical Solution of 3-Dimentional Poisson's Equation by Dimension Reduction Method(DRM), where the boundary conditions are Neumann.) です。

　以下,自分の過去論文をそのまま写します。なお,注は今回新たに書き入れたものです。　

                        Abstract

Up to the present time,iterative methods such as successive over-relaxation methods have been dominantly used to solve Poisson's equation numerically.

In many cases, however, it requires too much time to solve such partial differential equations by the iterative methods,even when a large-sized digital computer is used for calculation.

In this report,I describe a direct solving method of Poisson's equation in a rectangular parallelopiped domain where the gradient values are given on the boundary faces.

This method provides a precise solution and it is very economical of our calculation time.

１．序論

　本報告では,直方体領域におけるPoisson方程式のNeumann問題を,緩和法によらず直接的に解く数値解法を記述します。

　これには,1969年に,小倉¹⁾によって２次元のDirichlet問題の次元控除法による直接解法が与えられており,今回はそれを拡張して２次元,及び３次元のノイマン問題を解くことにします。

　なお,Dirichlet問題も同様に３次元に拡張することができます。

　この直接解法は,次元控除法(DRM法)と呼ばれ,Poisson方程式に特有の三重対角行列の固有値問題を解いて,固有ベクトルによる直交変換によって行列を対角化し,２次元,３次元の差分方程式を最終的に１次元の三重対角行列の係数の連立方程式に帰着させるものです。

　Poisson方程式のNeumann問題は,流体の速度が圧力のような変量の勾配で表わされ,境界の流速が与えられた場合に,その変量(例えば圧力)を導出する問題等において出現しますが,こうした直接解法が解の精度の向上や演算時間の短縮に非常に有効であると考えられます。

※  (注１)三重対角行列係数の連立方程式の解法と1次元Poisson方程式

   ｘ＝^ｔ(ｘ₁,ｘ₂,...,ｘ_n)を未知数とし,三重対角行列を係数とした

　ｎ元連立一次方程式とは,

　ａ₁ｘ₁＋ｃ₁ｘ₂＝ｙ₁,

　ｂ_kｘ_k-1＋ａ_kｘ_k＋ｃ_kｘ_k+1＝ｙ_k　(ｋ＝2,3,..,ｎ－1),

　ｂ_nｘ_n-1＋ａ_nｘ_n＝ｙ_n

なる方程式系で与えられる連立方程式のことをいいます。

これを行列で表現すると,定数項をｙ＝^ｔ(ｙ₁,ｙ₂,．．,ｙ_n)とし,

ｎ次正方行列Ａを,第１行目が(ａ₁,ｃ₁,0,...,0),第ｋ行目

(ｋ＝2,3,...,ｎ－1)が(0,..,ｂ_k,ａ_k,ｃ_k,0,..,0),第ｎ行目

が(0,...,0,ｂ_n,ａ_n)のように対角成分とその両側のみがゼロ

でない成分の行列＝三重対角行列であるとすれば,　

ｙ＝Ａｘと書けます。

　　

　これを解くのは簡単で,解法は大抵の数値計算の本に載っています。

ｄ₁＝ａ₁,ｅ₁＝ｙ₁として,ｄ_k+1＝ａ_k+1－ｃ_kｂ_k+1/ｄ_k,,

ｅ_k+1＝ｙ_k+1－ｅ_kｂ_k+1/ｄ_k(ｋ＝1,2,...,ｎ－1)と順次計算して,

その後ｘ_n＝ｅ_n/ｄ_n,ｘ_n-k＝(ｅ_n-k－ｃ_n-kｘ_n-k+1)/ｄ_n-kとする,

というアルゴリズムを用いれば解けます。

つまり,Ａｘ＝ｙと行列形式で書いたとき,Ａが三重対角行列なら

Ａ＝ＱＲ(Ｒは上三角行列),Ｑは直交行列)とＡをＱＲ分解して,

Ｒｘ＝ｅ＝Ｑ^-1ｙ＝^ｔ(ｅ₁,ｅ₂,...,ｅ_n)とします。

Ｒをｋ行目(ｋ＝1,2,...,ｎ－1)が(0,..,ｄ_k,ｃ_k,0,..,0),

第ｎ行目が(0,...,0,ｄ_n)(ただしｄ₁＝１とする)の上三角行列

とすると,

Ａ＝ＱＲ,かつｙ＝Ｑｅより,Ｑは第１行目が(1,0,...,0),

第ｋ行目(ｋ＝2,3,..,ｎ)が(0,..,ｂ_k/ｄ_k－1,1,0,..,0)で

与えられる下三角行列になるからです。

そして１次元のPoisson方程式は,φ(ｘ)を未知関数,ρ(ｘ)を与えられた周知の関数とすると,ｄ²φ/ｄｘ²＝ρ(ｘ)で与えられます。

これを差分方程式化すると,ｈをｘの差分として,

(φ_k+1－2φ_k＋φ_k-1)/ｈ²＝ρ_k　(ｋ＝1,2,...,ｎ)となります。

ただし,ｋ＝1のときはφ₀を,ｋ＝ｎのときはφ_n+1を境界条件によって

与える必要があります。

このため,直接境界値を与えたり(Dirichlet問題),境界での勾配:

ｆ₁＝(φ₁－φ₀)/ｈ,ｆ₂＝(φ_n+1－φ_n)/ｈを与える境界条件

(Neumann問題),あるいはこれらを混合した境界条件からφ₀と

φ_n+1を決めます。

例えば,Neumannン問題なら,これによって,

[(φ₂－φ₁)/ｈ－ｆ₁]/ｈ＝ρ₁より,(φ₂－φ₁)/ｈ²＝ρ₁＋ｆ₁/ｈ

が得られます。

そして,(－φ_n＋φ_n-1)/ｈ²＝ρ_n－ｆ₂/ｈです。

ρ₁をρ₁＋ｆ₁/ｈで,ρ_nをρ_n－ｆ₂/ｈで置き換えた列ベクトルを

ρと定義すれば,方程式はＡφ＝ρと行列形式で書けます。

このときｎ次の正方行列ｈ²Ａは第１行目が(－1,1,0,．． ,0),

第ｋ行目(ｋ＝2,3,...,ｎ－1)が(0,..,1,－2,1,0,..,0),

第ｎ行目が(0,..,0,1,－1)の三重対角行列になります。

ところが,実はこのＡ,あるいはｈ²Ａは固有値の１つがゼロの行列,

つまり非正則行列なので,逆行列が存在しませんから,単純に上記

の方法では解けません。

ただし,解ベクトルの1つの成分を適当に与えれば,漸化式的に他の

全ての成分を簡単に求めることができます。

それ故,これがNeumann問題でのPoisson方程式の解が微分方程式と

しても差分分方程式としても,なお定数だけの任意性を持つ所以と

なっています。

一応,三重対角行列を解くための,私の拙いFortranのソースコードを

下に入れておきます。

SUBROUTINE TRID(N,A,B,C,Y,X)

      REAL  A(N),B(N),C(N),Y(N),X(N)

      REAL  D(1000),E(1000)

      NN=N-1

      D(1)=A(1)

      E(1)=Y(1)

      DO 100 I=1,NN

      D(I+1)=A(I+1)-C(I)*B(I+1)/D(I)

      E(I+1)=Y(I+1)-E(I)*B(I+1)/D(I)

  100 CONTINUE

      X(N)=E(N)/D(N)

      DO 200 J=1,NN

      I=N-J

      X(I)=(E(I)-C(I)*X(I+1))/D(I)

  200 CONTINUE

      RETURN

　　　END 　　　　　　　　　　　　　　　　(注１終) ※

２．２次元ポPoisson方程式の解法(Neumann問題)

　矩形領域での２次元Poisson方程式は(∂²φ/∂ｘ²)＋(∂²φ/∂ｙ²)＝ρ(ｘ,ｙ) (ａ≦ｘ≦ｂ,ｃ≦ｙ≦ｄ)　...(2-1)で与えられます。

　Neumann問題での境界条件はｆ⁽¹⁾,ｆ⁽²⁾,ｇ⁽¹⁾,ｇ⁽²⁾を既知関数として

　ｘ＝ａで∂φ/∂ｘ＝ｆ⁽¹⁾(ｙ),ｘ＝ｂで∂φ/∂ｘ＝ｆ⁽²⁾(ｙ),

　ｙ＝ｃで∂φ/∂ｙ＝ｇ⁽¹⁾(ｘ),ｙ＝ｄで∂φ/∂ｙ＝ｇ⁽²⁾(ｘ)

...(2-2)で与えられるとします。

　まず,２次元のメッシュの個数がＮ_x×Ｎ_yであるとして差分ｈ_x,ｈ_y

を与えます：ｘ₀＝ａ,ｘ_Nx＝ｂ；ｘ_i－ｘ_i-1＝(ｂ－ａ)/Ｎ_x＝ｈ_x

　(ｉ＝1,2,..,Ｎ_x),ｙ₀＝ｃ,ｙ_Ny＝ｄ；ｙ_j－ｙ_j-1＝(ｄ－ｃ)/Ｎ_y＝ｈ_y

　(ｊ＝1,2,..,Ｎ_y) です。

　メッシュ(ｉ,j)の中心でφ,ρの値φ_i,j,ρ_i,jが定義されるとして

(2-1)式を差分化すると,

　[(φ_i+1,j－2φ_i,j＋φ_i-1,j)/ｈ_x²]

＋[(φ_i,j+1－2φ_i,j＋φ_i,j-1)/ｈ_y²]＝ρ_i,j

(ｉ＝1,2,...,Ｎ_x;ｊ＝1,2,...,Ｎ_y)　...(2-3)となります。

ただし,ｉ＝1のときにはφ_0,j,ｊ＝1のときにはφ_i,0,ｉ＝Ｎ_xのとき

にはφ_Nx+1,j,ｊ＝Ｎ_ｙのときにはφ_i,Ny+1の境界値が与えられなけれ

ばなりません。

これらは境界条件(2-2)で,ｆ⁽¹⁾(ｙ)をｆ⁽¹⁾_j,ｆ⁽²⁾(ｙ)をｆ⁽²⁾_j,ｇ⁽¹⁾(ｘ)をｇ⁽¹⁾_i,ｇ⁽²⁾(ｘ)をｇ⁽²⁾_iと離散化し,以下のように書き直すことで与えられます：

φ_0,j＝φ_1,j－ｈ_xｆ⁽¹⁾_j,φ_Nx+1,j＝φ_Nx,j＋ｈ_ｘｆ⁽²⁾_j,

φ_i,0＝φ_i,1－ｈ_yｇ⁽¹⁾_i,

φ_i,Ny+1＝φ_i,Ny＋ｈ_yｇ⁽²⁾_i

(ｉ＝1,2,..,Ｎ_x;ｊ＝1,2,..,Ｎ_y)...(2-4)

(2-4)を(2-3)に代入して,行列と列ベクトルによる連立方程式系に

まとめると,

φ_j+1－2φ_j＋φ_j-1＋(ｈ_y/ｈ_x)²Ａφ_j＝ρ'_j

(ｊ＝2,3,..,Ｎ_y－1) ...(2-5a),

φ₂－φ₁＋(ｈ_y/ｈ_x)²Ａφ₁＝ρ'₁  ...(2-5b),

－φ_Ny＋φ_Ny-1＋(ｈ_y/ｈ_x)²Ａφ_Ny=ρ'_Ny ...(2-5c)

となります。

ここに,φ_j＝^t(φ_1,j,φ_2,j,..,φ_Nx,j) (ｊ＝1,2,..,Ｎ_y),

ρ'_j＝^t(ρ_1,j＋ｆ⁽¹⁾_j/ｈ_x,ρ_2,j,..,ρ_Nx,j－ｆ⁽²⁾_j/ｈ_x)

(ｊ＝2,3,..,Ｎ_y－1),

ρ'₁＝^t(ρ_1,1＋ｆ⁽¹⁾₁/ｈ_x＋ｇ⁽¹⁾₁/ｈ_y ,ρ_2,1＋ｇ⁽¹⁾₂/ｈ_y,..,ρ_Nx,j

－ｆ⁽²⁾_j/ｈ_x＋ｇ⁽¹⁾_Nx/ｈ_y),

ρ'_Ny＝^t(ρ_1,Ny＋ｆ⁽¹⁾_Ny /ｈ_x－ｇ⁽²⁾₁/ｈ_y,ρ_2,Ny－ｇ⁽²⁾₂/ｈ_y..,ρ_Nx,Ny－ｆ⁽²⁾_Ny/ｈ_x－ｇ⁽²⁾_Nx/ｈ_y) 　 ...(2-6)です。

また,Ａは第１行目が(－1,1,0,..,0),第ｋ行目(ｋ＝2,3,..,Ｎ_x－1)

が(0,..,1,－2,1,0,..,0),第Ｎ_x行目が(0,..,0,1,－1)で与えられる

Ｎ_x×Ｎ_xの三重対角行列になります。...(2-7)

Ａの固有値方程式|λＥ－Ａ|＝0(Ｅは単位行列)は解析的に解くことができて,その固有値λ₁,λ₂,..,λ_Ｎxは,

λm＝－4sin²{(ｍ－1)π/(2Ｎ_x)}  ...(2-8)　で与えられます。

※  (注２):以下,三重対角行列Ａの固有値方程式|λＥ－Ａ|＝0 を具体的に解いてみます。

　簡単のために,Ａをｎ×ｎ行列とします。

固有値方程式|λＥ－Ａ|＝0 はλについてｎ次の代数方程式なので,正攻法で解くのはかなり困難です。

そこで,固有値方程式Ａｘ＝λｘを書き下してみます。

ｘ_k+1－(λ＋2)ｘ_k＋ｘ_k-1＝0 (ｋ＝2,3,..,ｎ－1)ですが,

通常の差分方程式(関数方程式＝漸化式)を解く手法に従って,

ｘ_k＝ξ^kとおいてこれを代入すれば,特性方程式:

ξ²－(λ＋2)ξ＋1＝0 が得られます。

これは,ξ＋1/ξ＝λ＋2と書けますから,ξ＝exp(iθ)とおけば,

2cosθ＝λ＋2,すなわち固有値λはλ＝2(cosθ－1)＝－4sin²(θ/2)

と書くことができます。

そして,ξ＝exp(iθ)がξ＋1/ξ＝λ＋2の解なら,ξ^-1＝exp(－iθ)もそうですから,ｘ_k+1－(λ＋2)ｘ_k＋ｘ_k-1＝0 の一般解はａ,ｂを任意定数としてｘ_k＝ａexp(iｋθ)＋ｂexp(－iｋθ)と表わすことができます。

ここで,この一般解が境界条件;ｘ₂－ｘ₁＝λｘ₁,

λ＝exp(iθ)＋exp(－iθ)－2 に適合する条件を求めると,　

ｂ＝ａexp(－iθ)となりますから,特にａ＝1,ｂ＝exp(iθ)

として,ｘ_k＝exp(iｋθ)＋exp{－i(ｋ－1)θ}とします。

これを－ｘ_n＋ｘ_n-1＝λｘ_nに代入すると,sin(ｎθ)＝0 を得ます。

したがって,sin(ｎθ)＝0 を満たす異なるｎ個のθの組を,

θ＝θ_ｍ＝(ｍ－1)π/ｎ (ｍ＝1,2,．．,ｎ)とすることができ

ます。

これらに対応して得られるｎ個の異なる固有値:

λ＝λ_ｍ＝2(cosθ_ｍ－1)＝－4sin²(θ_ｍ/2)が,

固有値方程式|λＥ－Ａ|＝0 のｎ個の根を示していると

考えられます。　(注２終わり)※

　Ａは実対称行列で固有値はすべて異なる値なので,直交行列Ｔ_Ｎx

　(すなわち,^tＴ_ＮxＴ_Ｎx＝Ｔ_Ｎx^tＴ_Ｎx＝Ｅ)が存在して,

　ＡはΛ_Ｎx＝^tＴ_ＮxＡＴ_Ｎx (Λ_Ｎxは対角成分がλ₁,λ₂,..,λ_Ｎxの

　対角行列)　...(2-9)と対角化できます。

　(※行列の下添字Ｎ_xは正方行列の次数を表わします。)

　Ｔ_Ｎxの陽な形は,後に付録(Appendix)で証明しますが,

　Ｔ_Ｎx＝(2/Ｎ_x)^1/2(ｔ₁,ｔ₂,ｔ₃..,ｔ_Nx);

　ｔ₁＝^t(1/√2,1/√2,1/√2,..,1/√2),

　ｔ₂＝(cos{π/(2Ｎ_x)},cos{3π/(2Ｎ_x)},cos{5π/(2Ｎ_x)},..

　,cos{(2Ｎ_x－1)π/(2Ｎ_x)}),

　ｔ₃＝(cos{2π/(2Ｎ_x)},cos{6π/(2Ｎ_x)},cos{10π/(2Ｎ_x)},..,

　cos{2(2Ｎ_x－1)π/(2Ｎ_x)}),..,

　ｔ_Nx＝(cos{(Ｎ_x－1)π/(2Ｎ_x)},cos{3(Ｎ_x－1)π/(2Ｎ_x)},

　cos{5(Ｎ_x－1)π/(2Ｎ_x)},..,cos{(2Ｎ_x－1)(Ｎ_x－1)π/(2Ｎ_x)})

　...(2-10a),

　すなわち,

　(Ｔ_Ｎx)_ｌｍ＝(1/Ｎ_x)^1/2　(ｍ＝1のとき),

　(Ｔ_Ｎx)_ｌｍ＝(2/Ｎ_x)^1/2cos{(2ｌ－1)(ｍ－1)π/(2Ｎ_x)}　

　(ｍ＝2,3,..,Ｎ_xのとき) ...(2-10b)

　です。

　そこで,Ｐ_j≡^tＴ_Ｎxφ_j,Ｒ_j≡^tＴ_Ｎxρ'_j (j=1,2,..,Ｎ_y)　...(2-11)

とおけば,(2-5)式は,

Ｐ_j+1－2Ｐ_j＋Ｐ_j-1＋(ｈ_y/ｈ_x)²Λ_ＮxＰ_j＝Ｒ_j(ｊ＝2,3,..,Ｎ_y－1)

　...(2-12a),

　Ｐ₂－Ｐ₁＋(ｈ_y/ｈ_x)²Λ_ＮxＰ₁＝Ｒ₁　...(2-12b),

　－Ｐ_Ny＋Ｐ_Ny-1＋(ｈ_y/ｈ_x)²Λ_ＮxＰ_Ny=Ｒ_Ny...(2-12c),

　すなわち,Ｐ_i,j+1＋{(ｈ_y/ｈ_x)²λ_ｉ－2}Ｐ_i,j＋Ｐ_i,j-1＝Ｒ_i,j

(ｊ＝2,3,..,Ｎ_y－1),

　Ｐ_i,2＋{(ｈ_y/ｈ_x)²λ_ｉ－1}Ｐ_i,1＝Ｒ_i,1,

　{(ｈ_y/ｈ_x)²λ_ｉ－1}Ｐ_i,Ny＋Ｐ_i,Ny-1＝Ｒ_i,1Ny　...(2-13)

　となります。

　これらはｉ＝1,2,...,Ｎ_xの各々を固定すると,ｊ＝1,2,...,Ｎ_y

に対して三重対角行列の係数を持った１次元の連立１次方程式系

です。そしてこれは次のように書けます。

すなわち,直交変換されたPoisson方程式の左辺の未知関数縦ベクトル

を,Ｐ'_i≡^t(Ｐ_i,1,Ｐ_i,2,..,Ｐ_i,Ny),右辺の既知関数縦ベクトル

をＲ'_i≡^t(Ｒ_i,1,Ｒ_i,2,..,Ｒ_i,Ny)とおき,行列Ａ'_iを

第１行目が({(ｈ_y/ｈ_x)²λ_ｉ－1},1,0,..,0),

第ｋ行目(ｋ＝2,3,..,Ｎ_x－1)が,

(0,..,1,{(ｈ_y/ｈ_x)²λ_ｉ－2},1,0,..,0),

第Ｎ_x行目が(0,,．,0,1,{(ｈ_y/ｈ_x)²λ_ｉ－1})

で与えられるＮ_y×Ｎ_yの三重対角行列とすれば,これを係数とし

未知数がＰ'_iの行列方程式; Ａ'_iＰ'_i＝Ｒ'_i (ｉ＝1,2,..,Ｎ_x)　

...(2-14) が得られます。

λ_i≠0 すなわち,i≠1のときには,左辺の行列Ａ'_iは正則であって,

逆行列が存在するので通常の１次元の三重対角行列係数の連立方程式

の解法を用いて解くことができます。

一方,λ_i＝0 すなわち,i＝1 のときには,左辺の行列Ａ'_iはＮ_y×Ｎ_y行列としてＡと同じ形になるので,非正則であるからＰ_1,1あるいはＰ_1,Nyを与えることによって,漸次Ｐ_1,jを計算できます。

こういう非正則の問題が生じるのは,Neumann問題を扱っているためであり解が定数だけの任意性を持つことに依ります。

　以上の過程で求めたＰ_jに対してφ_j≡Ｔ_ＮxＰ_jによってφ_jを求めれば

　２次元Poisson方程式のNeuman問題は解決されます。

　今日はここまでにします。

　以下は３次元ポPoison方程式 (Neumann問題) への拡張という論題

　へと続きます。

　この３次元の問題こそが,"成分が行列であるような行列＝直積

　＝テンソル積？"を導入するという私のオリジナルな発想を記述

　したものとなっています。

(※ 上記２次元問題は,元ネタの小倉さんの論文の紹介です。)

　参考文献：1)Masako Ogura："A direct solution of Poisson's equation by dimention reduction method." Journal of the Meteorogical Soc.Japan. Vol.47. No.4.(1969),pp319-322(日本気象学会誌、第47巻4号,1969) 2)矢嶋信夫、野木達男 著「発展方程式の数値解法」(岩波書店) (1977)　他

折り返しノイズ

前記事の続きです。

前記事ではｆ(ｔ)のFourier変換:Ｆ(ω)≡Ｆ(ｆ)において,|ω|＞ω₁のときはＦ(ω)＝0 であるという理想的なケースを考えました。

この場合には,Ｔ＜(π/ω₁)の間隔でサンプリングすれば完全に元の信号:f(ｔ)を復元できると書きました。

また,その式ではｆ(ｎＴ)のｎは－∞から∞の範囲となっていますが,実際には有限個しかサンプリングできないので,Shannonの公式はあくまで近似式であり,その分の誤差も生じることを書き忘れていました。。

そして,また実際にサンプリングするソースデータも理想的ではないケース,つまり,Ｔ＜(π/ω₁)の間隔Ｔでサンプリングしているのにも関わらず,元の音のソースは|ω|＞ω₁の領域にもゼロではないＦ(ω)の成分をも含むのが普通ですね。

理想的なケースなら,Ｆ(ｆ_T)＝(2π/Ｔ){∑_-∞^∞Ｆ(ω－2ｎπ/Ｔ)}の１つ１つのＦ(ω)の信号Ｆ(ω－2ｎπ/Ｔ)でのゼロでない帯域,－ω₁＜ω－2ｎπ/Ｔ＜ω₁は決して重なることはないです。

しかし,今の場合は|ω|＞ω₁の領域にゼロでないＦ(ω)の成分があるので,－ω₁＜ω－2ｎπ/Ｔ＜ω₁の部分に重なりができる結果,復元した信号の中で|ω|＞ω₁の高周波成分が,あたかもω₁＜ω－2ｎπ/Ｔ＜ω₁の成分であるかのように,ノイズとして入ってきます。

これをエイリアシング(Aliasing)とか,折り返しノイズと呼びます。

そして,こうしたノイズの発生を回避するためには,

１つには理想的ケースとするために,サンプリングの際に低周波のみをパス(スルー)し高周波をカットするフィルターであるローパスフィルター(Low-pass-filter)を通すという方法があります。

しかし,一般に完全なローパスフィルターというものはなく,これを実現するのは結構むずかしいことです。

そこで,別の手段として,サンプリング間隔Ｔをさらに小さくするという"オーバーサンプリング(Oversampling)"という手法を併用することが多いようです。

それでも,折り返しノイズを完全に除去するのは不可能で,ノイズを極限まで減らすことも技術的にはかなり困難なようです。

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

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

シャノンのサンプリング定理

　昨今はアナログ機器のデジタル化が進んでいますが,いくら細かくデジタル化しても結局はアナログを完全に再現することはできないのではないかと思います。

　集合の濃度という意味では,いくら稠密な集合であっても高々可算個のデータで連続無限を表現することは不可能ではないかと感じることはよくあります。

　この感覚とは別に"周波数の大きさに上限があるなら,時間的に連続な信号も上限の２倍以上の個数のサンプルがあれば再現(復元)可能である。"というShannon(シャノン)のサンプリング定理があります。

　以下では,それを説明してみようと思います。

　まず,準備として一般に関数ｆ(ｘ)のFourier変換:Ｆ(ｆ)を,

Ｆ(ｆ)≡Ｆ(ω)≡∫_-∞^∞f(ｘ)exp(－ｉωｘ)ｄｘで定義します。

　このとき,逆変換は,f(ｘ)＝[1/(2π)]∫_-∞^∞Ｆ(ω) exp(ｉωｘ)ｄωで与えられます。

また,周期がＴの周期的デルタ関数:δ_T(ｘ)を,

δ_T(ｘ)≡∑_-∞^∞δ(ｘ－ｎＴ)で定義します。

このとき,δ_T(ｘ)のFourier変換:Ｆ(δ_T)がやはり周期的デルタ関数になることがわかります。

すなわち,δ_T(ｘ)は周期Ｔの周期関数なので,これをFourier級数に展開すると,δ_T(ｘ)＝(1/Ｔ)∑_-∞^∞exp[i(2ｎπｘ/Ｔ)]となります。

そこで,Ｆ(δ_T)＝(1/Ｔ)∑_-∞^∞∫_-∞^∞[exp{－iｘ(ω－2ｎπ/Ｔ)}]ｄｘ＝(2π/Ｔ)∑_-∞^∞δ(ω－2ｎπ/Ｔ)と書けるからですね。

ｆ(ｔ)を帯域を制限された信号であるとします。

すなわち,Ｆ(ω)≡Ｆ(ｆ)において,|ω|＞ω₁のときはＦ(ω)＝0

であるとします。

例えば,音楽などでは人間の耳には,"周波数ｆ＝ω/(2π)が20kHz前後より大きい音＝非常に高い音(超音波と呼ばれ犬など他の動物には聞こえる音もあります)"は聞こえません7。

そこで,実質的にそれより大きい周波数部分は必要ないと見なされて,ＣＤでは22.05kHzより大きい周波数の信号はカットされています。

ここで,ある時間的周期Ｔごとにこの信号ｆ(ｔ)をサンプリングしたと考えると,サンプリングした信号をｆ(ｎＴ)として,全てのサンプル信号はｆ_T(ｔ)≡∑_-∞^∞ｆ(ｎＴ)δ(ｔ－ｎＴ)＝ｆ(ｔ)δ_T(ｔ)と表現できます。

これをFourier変換すると,Ｆ(ｆ_T(ｔ))＝Ｆ[ｆδ_T]

＝Ｆ[ｆ(ｔ)]*Ｆ[δ_T(ｔ)]となります。

最右辺の'*'は"たたみこみ(convolution)"を示しています。

そこで,Ｆ(ｆ_T)＝(2π/Ｔ)Ｆ(ω)*{∑_-∞^∞δ(ω－2ｎπ/Ｔ)}

＝(2π/Ｔ){∑_-∞^∞Ｆ(ω－2ｎπ/Ｔ)}となります。

もしも,ｆ_T(ｔ)のFourier変換Ｆ(ｆ_T)の周期(2π/Ｔ)が,周波数の幅である2ω₁よりも大きい:(2π/Ｔ)＞2ω₁ or Ｔ＜(π/ω₁)なら,

１つ１つの帯域幅を持つ信号:Ｆ(ω－2ｎπ/Ｔ)のゼロでない帯域:－ω₁＜(ω－2ｎπ/Ｔ)＜ω₁は決して重なることがないわけです。

そこで,元の信号の全体であるＦ(ω)が周期的に現われます。

ｆ(ｔ)は１周期分のFourier変換から逆変換により完全に決定できます。

具体的には,Ｆ(ｆ_T)は周期(2π/Ｔ)の周期関数なのでFourier級数に展開できて,Ｆ(ｆ_T)＝(2π/Ｔ){∑_-∞^∞Ｆ(ω－2ｎπ/Ｔ)}＝∑_-∞^∞ｃ_n exp(ｉｎＴω);ｃ_n ＝∫_-π/T^π/TＦ(ω)exp(－ｉｎＴω)ｄω＝Ｔｆ(－ｎＴ)です。

なぜなら,Ｆ(ω)は積分区間の外ではゼロなので積分区間を(-∞,∞)に拡張しても同じだからです。

そこで,Ｆ(ｆ_T)＝∑_-∞^∞Ｔｆ(ｎＴ)exp(－ｉｎＴω)であり,左辺は －π/Ｔ≦ω≦π/ＴではＦ(ω)に等しいので,

結局のところ,帯域内ではＦ(ω)＝∑_-∞^∞Ｔｆ(ｎＴ)exp(－ｉｎＴω)と展開されます。

これの逆変換により,元の連続信号ｆ(ｔ)が具体的に,

ｆ(ｔ)＝∑_-∞^∞[ｆ(ｎＴ)sin{π(ｔ/Ｔ－ｎ)}/{π(ｔ/Ｔ－ｎ)}]

という形で完全に得られます。

かくして,全信号は原理的にはサンプル信号ｆ(ｎＴ)だけで復元できることになります。

これを,"Shannonのサンプリング定理"と言います。

したがって,例えばＣＤの帯域が22.05 kHzだとすると,サンプリングには44.1kHzのデータがあれば十分ということになりますね。

つまり,帯域制限があるときには何らかの補間法に頼ることなく各時刻の各周波数ごとの振幅(大きさ)が計算によって復元できます。

つまり,どんな楽器のどんな音色やそれらの合奏であろうと,音の帯域周波数の２倍のデータをある場所でサンプルしておけば,その場所での全楽器からの音を原理的には完全に再現できるというわけです。

ただし,これらを本当に忠実に行なうのは結構むずかしい問題があると思われます。

これらの精度は,実際に復元演算を行なうデジタルコンピュータのチップを含むであろうＤＡコンバータ,録音機器,再生機器の性能にかかっており,

そうしたことを行なうプロセスでは,必然的にノイズやエラーを伴なうからです。

また,実はカットした帯域外の人間の耳には聞こえないはずの音が弾性振動の波として違った作用を人体に及ぼす結果として,人間の快,不快の感覚に関与しているかも知れません。

参考文献；松下泰雄 著「フーリエ解析」(培風館)

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

人気blogランキングへ　← クリックして投票してください。(１クリック＝１投票です。１人１日１投票しかできません。）

← クリックして投票してください。(ブログ村科学ブログランキング）

    物理学

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

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

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

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

　例えば容積がＶ(ｍ³)の部屋の中央に石油ストーブなどがあって,それからの一酸化炭素(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(ブログ・ランキングの投票）↑ここをクリックすると投票したことになります。

非線形最小二乗法（JEA式の作成過程)

　今日は窒素酸化物濃度などを計算する解析的な低煙源拡散式のパラメータの推定のために用いたことがある非線形最小二乗法について解説してみたいと思います。(直角風のJEA式)

　地上の道路上を走行する自動車からの窒素酸化物(NOx)などの汚染物質排出のソースを,有効煙源高さ(effective height):ｈe＝0 (ｍ)のｙ軸上で単位長さ当たりの定常的な排出強度がＱ(Ｎｍ³/(ｍｓ)の無限線煙源としてモデル化します。

　そして,道路に直角のｘ方向に風が吹いている設定で,風速も拡散係数も高さｚ(ｍ)のベキ乗則に従って鉛直上方に向かって増加するというモデルを想定します。

　このときの拡散方程式の解である汚染物質濃度を,垂直距離ｘ(ｍ)と高さｚ(ｍ)の関数と考えてＣ(ｘ,ｚ)で表わすと,これはＡをＱに比例するあるパラメータとして,Ｃ(ｘ,ｚ)＝Ａｘ^-sexp(－Ｂｚ^p/ｘ)なる形の式(Robertsの式)で表わされる,ことがわかります。

　様々な環境の下で数回に渡って行なわれた実際の実験時の煙源長さはもちろん有限なのですが,とりあえず,これを無限大長さで近似します。

　　　

　実はガウスの誤差関数でもって有限長さの効果を取り入れることもできるのですが,今日の話題では割愛します。

　実験はトレーサーガスを地上にある有限長さの直線状のパイプに均等間隔で開けた穴から拡散させるというものです。

そのパイプ状の発生源をｙ軸としたとき,ｎ個の観測点座標:(ｘ_i,ｚ_i)(ｉ＝1,2,..ｎ)において,風向が直角に近い環境のときの実測濃度:ｃ_iと先に設定した低煙源拡散式による計算値:Ｃ_i＝Ｃ(ｘ_i,ｚ_i)とを比較することによって,逆にその式のパラメータＡ,ｓ,Ｂ,ｐを推定します。

　具体的には,Ｃ(ｘ,ｚ)＝Ａｘ^-sexp(－Ｂｚ^p/ｘ)の右辺をｆ(Ａ,ｓ,Ｂ,ｐ；ｘ,ｚ)と書いて,Ｃ_i＝ｆ(Ａ,ｓ,Ｂ,ｐ；ｘ_i,ｚ_i)として,この計算値と実測値ｃ_iの誤差の二乗和:Ｓ(Ａ,ｓ,Ｂ,ｐ)≡∑(Ｃ_i－ｃ_i)²を最小にするパラメータＡ,ｓ,Ｂ,ｐを求めるわけです。

　これは,次に述べる非線形最小二乗法のパラメータ数が４個のときの例になっています。

　ここでは,より一般的に未知パラメータがｒ個あるとし,それらを順にＡ₁, Ａ₂,..,Ａ_rとします。

　上のケースではｒ＝4で,Ａ₁＝Ａ,Ａ₂＝ｓ,Ａ₃＝Ｂ,Ａ₄＝ｐです。

　そして,Ｃ_i＝Ｃ(ｘ_i,ｚ_i)＝ｆ(Ａ₁,Ａ₂,..,Ａ_r；ｘ_i,ｚ_i)とし,Ｓ(Ａ₁,Ａ₂,..,Ａ_r)≡∑(Ｃ_i－ｃ_i)²と定義するわけです。

　誤差の二乗和Ｓが最小となる条件は,

　通常の線形回帰の計算の場合と同様,必要条件として

　Ｆ_i(Ａ₁,Ａ₂,..,Ａ_r)≡∂Ｓ/∂Ａi＝2∑(Ｃ_k－ｃ_k)∂Ｃ_k/∂Ａi＝0 (ｉ＝1,2,..,ｒ)　で与えられます。

　そして,この r 個の連立方程式を解くことにより未知数Ａ₁,Ａ₂,..,Ａ_rを求めることが主目的となります。

　これらの方程式は一般に非線形ですから,こうした非線形回帰によるパラメータ推定の方法を非線形最小二乗法と呼ぶことにします。

　具体的には,Ｆ_i(Ａ₁,Ａ₂,．．Ａ_r)＝0 (ｉ＝1,2,．．ｒ)を線形近似することにより多変数ニュートン法(Newtonian method)を実行します。

　そこで,まず連立方程式を線形近似します。

　すなわち,初期値としてＡ_j＝Ａ_j⁰ を適当に与えた後に,

　0 ＝ Ｆ_i(Ａ₁,Ａ₂,..Ａ_r)

　 ＝Ｆ_i(Ａ₁⁰,Ａ₂⁰,..,Ａ_r⁰)＋∑(∂Ｆ_i/∂Ａ_j)|_Ａ＝Ａ0(Ａ_j－Ａ_j⁰)

　と近似します。

これを行列形式で書くため,Ｄという行列をＤ＝(ｄ_ij)≡(∂Ｆ_i/∂Ａ_j)で定義し,特にＤ₀＝(ｄ_ij⁰)≡(∂Ｆ_i/∂Ａ_j)|_Ａ＝Ａ0とします。

一方,列ベクトルＡをＡ≡^t(Ａ₁,Ａ₂,..,Ａ_r)で定義し,特に

Ａ₀≡^t(Ａ₁⁰,Ａ₂⁰,..,Ａ_r⁰)とします。

さらに,Ｆ_i⁰≡Ｆ_i(Ａ₁⁰,Ａ₂⁰,..,Ａ_r⁰)を成分とする同様な列ベクトルをＦ₀と定義します。Ｆ₀≡^t(Ｆ₁⁰,Ｆ₂⁰,..,Ｆ_r⁰)です。

こうすれば,先の線形近似は 0＝Ｆ₀＋Ｄ₀(Ａ－Ａ₀)と表わされます。

これを単純に解くと,Ａ＝Ａ₀－Ｄ₀^-1Ｆ₀と近似されることになります。

ここでＤ₀^-1はＤ₀の逆行列です。

得られたＡ＝Ａ₀－Ｄ₀^-1Ｆ₀を,改めてＡ₀として代入してＤ₀^-1Ｆ₀を計算し,Ａ＝Ａ₀－Ｄ₀^-1Ｆ₀が収束するまで繰り返します。

すなわち,Ａ_m+1＝Ａ_m－Ｄ_m^-1Ｆ_mの漸化式において誤差|Ａ_m+1－Ａ_m|の相対値が十分小さくなるまで繰り返し計算します。

具体的には,Ａ₁＝Ａ₀－Ｄ₀^-1Ｆ₀,Ａ₂＝Ａ₁－Ｄ₁^-1Ｆ₁,..,Ａ_m+1＝Ａ_m－Ｄ_m^-1Ｆ_mを用います。ただし,Ｄ_k＝(ｄ_ij^k)≡(∂Ｆ_i/∂Ａ_j)|_Ａ＝Ａkです。

ｍ→ ∞では,|Ａ_m+1－Ａ_m|→ 0,それ故Ｆ_m→ 0 となることが予想されるため,この反復計算を実行するわけです。

こうすれば,ｍ→ ∞ でのＡ_m＝^t(Ａ₁^m,Ａ₂^m,..,Ａ_r^m)の極限値としてＦ_i(Ａ₁,Ａ₂,..,Ａ_r)＝0 (ｉ＝1,2,．．．ｒ)を満たすＡ＝^t(Ａ₁,Ａ₂,..,Ａ_r)が得られるはずです。

さらに,実際には収束を速くするために加速係数:ｗを与えて,Ａ_m+1＝Ａ_m－ｗＤ_m^-1Ｆ_mとする方法などを用いた方がいいと思います。

(※ここではｗを加速係数であると称してはいますが,それは広い意味でありｗ＞1を意味するわけではありません。実際,計算ではｗとして2.0や逆に0.01など収束が悪い場合には,さまざまな値で試行錯誤しました。)

実際のｒ＝4 の計算では,かつてFortranを使ってプログラミングしましたが,初期値を適切に取ると各環境のケースについて大体十数回の繰り返し計算でパラメータの推定値が得られました。

そして濃度の実測値とその推定パラメータによる計算値との相関係数としては,0.9 前後の値が得られ,回帰係数の値も0.8 から 1.2程度となって,仮定した計算式が良い近似になることがわかりました。

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

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