311 .数値計算・調和解析・離散数学

2007年7月16日 (月)

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

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

 

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

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

 

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

すなわち,各枝は容量制限Lm≦xm≦Um (Lm,Umは±∞,正,負,0のいずれでも良い)を有するが,この他にi=-jm,..,0,..,kmの各値について,tm,iが存在してtm,i-1≦xm≦tm,iの各xmに対して,単位コストcm,iを持つ。なお,tm,00 である。

mに関わる総コストは,xm 0 のときは,tm,lm-1≦xm≦tm,lmならばΣk=1lm-1m,k(tm,k-tm,k-1)+cm,lm (xm-tm,lm-1)(3-1a)であり,xm 0 のときは,tm,-lm-1≦xm≦tm,-lmならばΣk=0lmm,-k(tm,-k-1-tm,-k)+cm,-lm (xm-tm,-lm-1)(3-1b)である。

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

 

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

基底枝em(i,j)に対しては,tm,k-1≦xm≦tm,kにおける単位コストcm,kについて双対基底として,Πj+cm,k=Πi(3-2)となるように、{Πn}(n=1,2,..,N)の組を計算する。

 次に,m(i,j)が非基底のとき,m=tm.k-1ならばDm=Πi-Πj-cm,k ,xm=tm.kならばDm=Πj-Πim,kとおく。 

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

 

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

 

つまり,全てのmがDm0 なら,コスト関数においてxmの±Δmの変化と共に変化し得る量Σcm,kmm(2-2)式と同様に変形されて,Δm>0 の変化に対して必ず非減少となるからである。

そこで次に,得られた暫定的最小コスト解において,非基底枝mに対してmが下界tm.km-1に等しい場合には,これを上界と見なして,新しいDmの値としてのDm*をDm*=Πj-Πim,km-1で与える。

 

m0 より-Dm=Πj-Πim,km0 であり,単位コストm,kはkに関して単調増加であるからDm*≦-Dmであるが,もしDm*0 ならmは基底に入る枝の候補となる。 

同様にmが上界tm.kmに等しい場合には,これを下界と見なして,新しいDmの値としてのDm*をDm*=Πi-Πjm,km+1で与える。

 

m0 より-Dm=Πi-Πjm,km0 であり,単位コストm,kはkに関して単調増加であるから,やはりDm*≦-Dmであるが,もしDm*0 ならmは基底に入る枝の候補となる。 

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

 

なぜなら,(3-1)式を,m=1からMまで加えた総コストはΣmk≦km-1m,k(tm,k-tm,k-1)+cm,kmm,km-1}+Σm∈Umm,km-Σm∈Lmm,km-1+ΣnΠnn-Σ¬em∈TmΔm(3-3)と書けるから,Dm≦0 によりΔm>0 に対しては既に最小であり,上述の操作による非基底枝の境界を越えるΔm<0 の方向へのmの変化に対しては,コスト関数を(3-3)式の形に変形した場合,xmの変化Δm*=-Δm0 についても,なお更新されたDmの値Dm*が全てゼロ以下となるからである。

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

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

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

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

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

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

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

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

 

流量は,全て非負の値にするため,新しい流量変数をxm'≡xm-Lm (m=1,2,..,M)(A-1)で与え,上界をUm-Lmとする。

(1-1)より,rn=Σem∈Anm-Σem∈Bnmであるから,新しい供給量rn'をrn'=Σem∈Anm'-Σem∈Bnm'=rn-Σem∈Anm+Σem∈Bnm (n=1,2,..,N)(A-2)と変換して,rN+1'=0 とおく。

 人為枝eM+iはri'>0 のとき(i,N+1)で,ri'≦0 のとき(N+1,i)とする。

 

また,全ての枝emの単位コストをpm(m=1,2,..,M,M+1,..,M+N)として,emが実枝ならpmに 0 を人為枝なら 1 を与えることにより1つの線型輸送問題ができる。

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

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

 条件:Σm∈Anm'-Σem∈Bnm'=rn'(n=1,2,..,N,N+1) (A-3); 0≦xm'≦ Um-Lm (m=1,2,..,M)(A-4a);0≦xm'(m=M+1,..,M+N)(A-4b)

 

 目的:z=Σm=1M+Nmm'(A-5)を最小化する。

 

 ただし,pm0 (m=1,2,..,M)(A-6a);m1(m=M+1,..,M+N)(A-6b)である。

補助問題の初期可能木解はN個の基底枝を人為枝eM+1,eM+2,..,eM+Nとして,それらの流量の初期値をri'>0 ならxM+i'= ri'(A-7a),ri'≦0 ならxM+i'=-ri'(A-7b)と設定し,M個の非基底枝em(m=1,2,..,M)については,流量をxm'=0 (m=1,2,..,M)(A-8)と設定すれば得られる。

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

以下,z=0 の場合のみを考える。

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

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

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

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

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

参考文献 

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(V.フバータル:阪田省二郎,藤野和雄,田口 東 訳)「線形計画法」(上)(下)(啓学出版)(1988)

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

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

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

 

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

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

 

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

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

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

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

 

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

 

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

 

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

 

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

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

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

 

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

 

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

 

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

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

 

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

 

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

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

 

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

  

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

 

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

 

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

 

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

 

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

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

 

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

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

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (0) | トラックバック (0)

2007年7月15日 (日)

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

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

 

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

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

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

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

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

節点iの先行点piとは,iから根までの単純経路上にあって,深さが(di1)の節点のことである。p1=0 と定義する。

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

 

一方,fを節点集合V={1,2,..,N}からVの上への1対1の写像であるとする。Ui≠φでf(j)=iのとき,{f(k):j+1≦k≦j+|Ui|}=Ui (ただし,|A|は有限集合Aの元の個数)が成立すれば,この写像fを1つの継続点配列と呼ぶ。

1つの継続点配列が与えられたとき,節点iの糸siとは,この配列でiに続く節点,つまりf(j)=iならf(j+1)のことである。配列の最後の糸は根1であると定義する。

iに根付く極大木の最後の節点をniとすると,niは次の条件のうちの1つを満足しなければならない。:(a)sni=1(b)psniは節点iから根までの道の上にある。

(注)つまり,sni≠1ならpsni≠0 であり,sniから根までの道の上でのsniの先行点がpsniなので,psniはもちろんsniから根までを結ぶ道の上にある。

 

 そしてsniはniから根1までをつなぐ|Uni|個の点のうちのni以外のどれか1つである。

 

 しかもniはiに根付く極大木の最後の節点だからsniはi→ni→1という道の上にあり,psniはsni→1なる道の上にあるから,結局psniはiと根1を結ぶ道の上にあることになる。

次にeklがiとpiを結ぶ枝であるとする。

 

kl(pi,i)のとき,mi=ki ,ekl=(i,pi)のとき,mi=-k と定義する。ekl上の流量はαiで表わすことにする。

基底に入る枝をeφ(u,v)とするとき,更新手順は以下の「アルゴリズム1」および「アルゴリズム2」に示す通りである。

(1)アルゴリズム1

 

① 初期化

. γ1とγ2を初期化する。

φ=Lφならγ1=-1,γ2=1;

φ=Uφならγ1=1,γ2=-1とおく。

b.Δを臨界値としておく。

Δ=Uφ-Lφ とおく。

1とδ2を初期化する。

δ1=u2=vとおく。

min(du,dv)より大きい深さラベルを持つuからvへの部分経路上の最大許容流量変動Δを決定する。

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

u-dv0 ならd=dv-du ,μ=2とおく。du-dv=0 なら③へ進む。du-dv>0 ならd=du-dv ,μ=1とおく。

b.δ,γと計数を設定する。

δ=δμ,γ=γμ,k=1とする。

.増加流か減少流か?

δγ>0 ならeへ進む。

.増加流

|mδ|-αδ>Δなら,δ=pδとしてfに進む。さもなければ,Δ=U|mδ|-αδ,λ=δ,δ=pδとしてfに進む。

.減少流

αδ-L|mδ|>Δなら,δ=pδとする。さもなければ,Δ=αδ-L|mδ|,λ=δ,δ=pδとする。

.δの深さ=min(du,dv)

k<dならk=k+1とおきcに進む。

. δμの再設定

δμ=δとおく。

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

.道の終わりか?

δ1=δ2なら,w=δ1として終了する。さもなければk=1とおく。

.δのロード

δ=δk とおく。

.増加流か減少流か?

δγk0 ならeへ進む。

.増加流

|mδ|-αδ>Δなら,δk=pδとしてfに進む。さもなければ,Δ=U|mδ|-αδ,μ=k,λ=δ,δk=pδとしてfに進む。

.減少流

αδ-L|mδ|>Δなら,δk=pδとする。さもなければ,Δ=αδ-L|mδ|,μ=k,λ=δ,δk=pδとする。

.両側を調べたか?

k=1ならk=2として,bへ進む。さもなければaへ進む。

「アルゴリズム1」の終了時には,最大許容流量変動Δ,出る枝e|mλ|が定まり,uからvへの道はuからwの道とvからwの道によって与えられる。

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

(2)アルゴリズム2

① 全更新が必要か否か?

Δ≠Uφ-Lφなら③へ進む。

②更新流の再設定

.閉路に対する初期設定

δ1=u2=v,k=1とする。

.枝eφの流量更新

φ=Uφならxφ=Lφとする。さもなければxφ=Uφとする。

.パラメータの初期化

δ=δk,γ=γkΔとおく。

d.節点かwか?

δ=wならhへ進む。

.増加流か減少流か?

δ0 ならβ=1,mδ<0 ならβ=-1とおく。

f.更新流の再設定

αδ=αδ-γβとおきかえる。

g.次の先行点

δ=pδとしdに進む。

h.k=1ならk=2とおきcに進む。さもなければ終了する。

③初期化

a.双対変数の変更

σ=Dφとおく。

b.入る枝の更新流量は?

φ=Uφなら,Δ'=Uφ-Δ,σ=-σとする。さもなければ,Δ'=Lφ+Δとする。

C.q,q'の初期化

μ=1ならq=u,q'=v,ε=-φ,σ=-σとする。さもなければq=v,q'=u,ε=φとする。

d.i,j,γの初期化

i=q,j=-p,γ=γkΔとおく。

e.λの先行点のセーブ

λ'=pλとしておく。

f.q'の深さラベルをセーブ

F=dq'1としておく。

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

.双対変数の更新

Πi=Πi+σとおく。

.流量のセーブ

α'=αiとしておく。

c.流量の更新

αi=Δ'とする。

d.流れの方向のセーブ

i0 ならβ=1, mi<0 ならβ=-1とおく。

e.枝のセーブ

m=|mi|と枝番号をセーブする。

 

.枝ラベルの更新

i=εとする。

g. iの深さdi ,およびFとdiの差をセーブする。

'=di,L=F-diとしておく。

h. iの深さの更新

i=Fとする。

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

a. iから始める。

z=i,x=siとおく。

b. xは1つの継続点か否か?

x≦F'なら⑥へ進む。

c.双対変数とxの深さの更新

Πx=Πx+σ,dx=dx+Lと更新する。

d.続く糸

z=x,x=sxとしてbへ進む。

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

. jから始める。

r=jとおく。

b. rは逆糸か?

riならば⑦へ進む。

.続く糸

r=srとしてbに進む。

⑦先行点順列と糸の更新

.最終更新か?

i=λならiに進む。

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

Δ'=α'-γβとおきかえる。

c. jの枝ラベルを決める。

ε=-βmとする。

d. rとzの糸の更新

r=x,sz=jとおく。

e. iをセーブする。

r=iとおく。

f. iの更新

i=jとおく。

g. jの更新

j=pjとおく。

h. iの先行点の更新

i=r,F=F+1として④へと進む。

i. rとzの糸の更新

r=x,sz=sq'とする。

j. q'の糸の更新

q'=qとする。

k. qの先行点の更新

q=q'とする。

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

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

a.閉路に対する設定

μ=1ならδ1=λ',δ2=q',k=1とおく。さもなければ,δ1=q',δ2=λ',k=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(V.フバータル:阪田省二郎,藤野和雄,田口 東 訳)「線形計画法」(上)(下)(啓学出版)(1988)

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

 

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

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

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (0) | トラックバック (0)

2007年7月13日 (金)

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

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

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

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

 

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

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

特にファインマン(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.     序文 

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

 

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

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

 

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

各枝em(i,j)は,節点iから節点jへの向きを持った流れを持つとする。(m=1,2,..,Mで,i,jは1,2,..,Nの選ばれた異なる2つの番号である。)

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

 

n>0 ならnはrnに等しい供給を与える供給点である。rn<0 ならnは-rnに等しい需要を与える需要点である。またrn=0 ならnは中間点である。つまり,rnは節点nにおける代数的供給量である。

任意の枝em(i,j)を通過する流れがあるとき,単位流量当たりのコストがcmであるとする。

 

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

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

各節点nに対し,nから出て行く枝の集合をAn,nに入ってくる枝の集合をBnとするとき,各節点n(n=1,2,..,N)に対して,Σem∈Anm=Σem∈Bnm+rn(1-1)が成立しなければならない。

 

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

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

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

 

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

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

 

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

2-1.木と可能木解

①節点uとvの間の道とは,節点w1,w2,..,wkと枝e1,e2,..,ek-1から成るネットワークで,u=w1,v=wk かつ枝eiの両端点がwi,wi+1であるようなものである。

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

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

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

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

 

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

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

 

 つまり、枝emに対しLmをその流量の下界,Umを流量の上界とするとき,emがTに属さないならxm=Lm,またはxm=Umであって,全てのnについて,条件(1-1)を満足する{xm}の組(m=1,2,..,M)をTに対する可能木解という。

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

 

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

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

 求めるアルゴリズムに到る最初のステップは,各節点nに対する双対変数 Πnとして全ての基底枝em(i,j)に対し,Πj+cm=Πi(2-1)となるような変数:{Πn}(n=1,2,..,N)の組を計算することである。

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

 次に,em(i,j)が非基底枝の場合(すなわちemがTに属さない場合:¬em∈T)には,xm=LmのときDm=Πi-Πj-cmとおき,これらの枝の集合をL,xm=UmのときDm=Πj-Πi+cmとおき,これらの枝の集合をUとする。

 

 ただしLm=Umのときには,どちらか一方を適当に選ぶ。

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

 なぜなら,Σm∈Anm=Σm∈Bnm+rn (1-1)の条件を満たしながら,非基底についてxm=Lm+Δm,またはxm=Um-Δmm>0)とxmを変化させたとき,(それによって基底のxmは自動的に変化する。)

Σm∈Tj-Πi+cm)xm=0 だから,総コストは,Σm=1Mmm=Σm=1Mj-Πi+cm)xmm+Σm=1Mi-Πj)xm=Σem∈Tj-Πi+cm)xm+Σ¬em∈Tj-Πi+cm)xm+ΣiΠiΣem∈Aim-ΣjΠjΣem∈Bjm=Σ¬em∈Tj-Πi+cm)xm+ΣnΠnem∈Aim-Σem∈Bjm)=Σem∈Umm-Σem∈Lmm+ΣnΠnn-Σ¬em∈TmΔm (2-2)

 そこで,全てのDmがDm0 のとき,任意の非基底em上のΔm>0 の変化に対して,Σcmmは現在の状態より必ず非減少となるため,現在の可能木解が1つの最小コスト解となる。

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

 

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

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

 

 基底に入る枝emの流量xmの変化を±Δとすると,Δ>0 が大きいほどコストの減少は大きいわけであるが,xmの±Δだけの変化に対して閉路内の基底枝の全てがそれぞれ±Δずつ変化するときに限って,(1-1)式のΣem∈Anm=Σem∈Bnm+rnが満足される。

閉路に属する全ての基底についての容量制限Lm≦xm≦Um(em∈T)を満足しつつ,最大の値となるようにΔの値を定めることができる。

 

それにつれて,xm=Lm,またはxm=Umとなって,基底から出て非基底となる1つの枝が得られる。

 

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

入る枝に対してDm0 ,出る枝に対してDm=0 であるから,このプロセスでコストは減少する。

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

 

Tに入る枝をeφ,Tから出る枝をeψとすると,新しい極大木はT+eφ-eψと書ける。

 

T+eφ-eψから枝eφを除いたネットワーク(T+eφ-eψ)-eφ=T-eψは,枝eφの端点を1つずつ含む2つの互いに素な木から成る。

φ(u,v)のとき,節点uを含む木をTu,vを含む木をTvとする。

 

φの単位コストがcφなら,σ=Πv-Πu+cφとおいて既知のΠnから更新されたΠnの値Πn*を次式によって求める。

 

Πn*=Πn+σ(n∈Tu),Πn*=Πn (n∈Tv)(2-3)である。

 

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

 

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

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

 

参考文献

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(V.フバータル;阪田省二郎,藤野和雄,田口 東 訳)「線形計画法」(上)(下)(啓学出版)(1988)

 

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

 

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

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

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (0) | トラックバック (0)

2007年5月21日 (月)

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

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

 

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

 

つまり,△x()=δ()となるGreen関数;D()を求める

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

します。

 

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

特に,球対称なGreen関数は,D()=-1/(4πr)(ただしr=||)で

与えられます。

 

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

 

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

そうして得られたD()を用いて,φ(x)=∫D()ρ()dとおけば,φは確かに△φ=∇2φ=ρの解になります。

 

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

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

 

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

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

 

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

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

 

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

 

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

ところです。 

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

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

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

 

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

 

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

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

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

 

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

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

 

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

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

 

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

  

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

  

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

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

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (0) | トラックバック (0)

2007年5月20日 (日)

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

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

 

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

(∂2φ/∂x2)+(∂2φ/∂y2)+(∂2φ/∂z2)=ρ(x,y,z)

(a≦x≦b,c≦y≦d,e≦z≦f)...(3-1)で与えられます。

 

Neumman的境界条件は,(∂φ/∂x)|x=a=f(1)(y,z),

(∂φ/∂x)|x=b=f(2)(y,z),(∂φ/∂y)|y=c=g(1)(x,z),

(∂φ/∂y)|y=d=g(2)(x,z),(∂φ/∂z)|z=e=h(1)(x,y),

(∂φ/∂z)|z=f=h(2)(x,y) ...(3-2) です。

 

 直方体を,各辺がhx,hy,hzであるようなNx×Ny×Nz個の微小セル

分割します。

 

  x0=a,xNx=b;xi-xi-1=(b-a)/Nx=hx(i=1,2,..,Nx);

  y0=c,yNy=d;yj-yj-1=(d-c)/Ny=hy(j=1,2,..,Ny);

  0=e,zNz=f;zk-zk-1=(f-e)/Nz=hz(k=1,2,..,Nz)

 

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

  

  [(φi+1,j,k-2φi,j,k+φi-1,j,k)/hx2]

 +[(φi,j+1,k-2φi,j,k+φi,j-1,k)/hy2]

 +[(φi,j,k+1-2φi,j,k+φi,j,k-1)/hz2]=ρi,j,k

 

 (i=1,2,...,Nx;j=1,2,..,Ny;k=1,2,..,Nz).

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

 

 φ0,j,kNx+1,j,ki,0,ki,Ny+1,ki,j,0i,j,Nz+1は境界条件

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

 ます。

 

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

変換すると,

 

Φk+1-2ΦkΦk-1+{(hz/hy)2A~+(hz/hx)2B~}ΦkΡ'k

(k=2,3,..,Nz-1)...(3-4a), 

Φ2Φ1+{(hz/hy)2A~+(hz/hx)2B~}Φ1Ρ'1 (3-4b),

ΦNzΦNz-1+{(hz/hy)2A~+(hz/hx)2B~}φNzΡ'Nz 

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

 

Φkt(φ1,k,φ2,k,..,φNy,k),

φj,kt1,j,k2,j,k,..,φNx,j,k),

Ρ'kt(ρ'1,k,ρ'2,k,..,ρ'Ny,k) 

と定義します。

   

ρ'j,kt1,j,k+f(1)j,k/hx+g(1)1,kδj1/hy-g(2)1,kδjNy/hy

+h(1)1,jδk1/hz-h(2)1,jδkNz/hz2,j,k

+g(1)2,kδj1/hy-g(2)2,kδjNy/hy

+h(1)2,jδk1/hz-h(2)2,jδkNz/hz,..,ρNx,j,k-f(2)j,k/hx

+g(1)Nx,kδj1/hy-g(2)Nx,kδjNy/hy+h(1)Nx,jδk1/hz

-h(2)Nx,jδkNz/hz)

 

(j=1,2,..,Ny;k=1,2,..,Nz)...(3-5)です。

 

 (Nx×Ny)次の正方行列A~,B~はNx次の単位行列Eと,

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

 表わすことができます。

 

 A~は第1行目が(-E,E,0,..,0),第k行目(k=2,3,..,Ny-1)

 が(0,..,E,-2E,E,0,..,0),第Nx行目が(0,..,0,E,-E)で

 与えられるNx次正方行列を成分とするNy×Ny

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

 また,B~は対角成分が全てNx次正方行列Aでそれ以外の成分は

 ゼロ のNx次正方行列を成分とするNy×Ny

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

  ただし,Nx次正方行列Aは次の通りです。

 A~を対角化する(Nx×Ny)次の直交行列をS~Nyとすると,

 tS~NyA~S~NyはEをNx次の単位行列として,対角成分が,

 σ1E,σ2E,..,σNyEの対角行列になります。 ...(3-8)

 ここに,

 σm=-4sin2{(m-1)π/(2Ny)} (m=1,2,..,Ny)...(3-9)

 

  (S~Ny)il,jm=(1/Ny)1/2δij (m=1;l=1,2,..,Ny),

  (S~Ny)il,jm=(2/Ny)1/2cos{(2l-1)(m-1)π/(2Ny)}δij

  (m=2,3...,,Ny;l=1,2,..,Ny) (i,j=1,2,..,Nx) .(3-10)

 

です。

 

 また,前記事の(2-9)によって,ΛNxtNxATNxNxは対角成分

 がλ12,..,λNxの対角行列) ...(3-11)であり, 

 λm=-4sin2{(m-1)π/(2Nx)} (m=1,2,..,Nx)です。

 

 T~Nyを対角成分が全てTNxのNy次対角行列 (3-12)と考えれば,

 tT~NyB~T~Nyは対角成分が全てΛNxのNy次対角行列

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

 

 (T~Ny)il,jm=(1/Nx)1/2δlm (j=1;i=1,2,..,Nx),

 (T~Ny)il,jm=(2/Nx)1/2cos{(2i-1)(j-1)π/(2Nx)}δlm

 (j=2,3,..,Nx;i=1,2,..,Nx) (l,m=1,2,..,Ny)

 ..(3-14)です。

 

 以上から,A~の直交変換:t(S~NyT~Ny)A~S~NyT~Nyは,対角成分が

 σ1E,σ2E,..,σNyEの対角行列になります。..(3-15),

 

 一方,B~の直交変換:t(S~NyT~Ny)B~S~NyT~Nyは,対角成分が全て

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

それ故,Q~≡S~NyT~Nyと置けば,Q~による直交変換によって,

A~,B~は同時に対角化できて,

 

tQ~{(hz/hy)2A~+(hz/hx)2B~}Q~は,対角成分が,

 [(hz/hy)2σmE+(hz/hx)2ΛNx]

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

 

 この得られた対角行列をΩとし,kt~Φk,ktQ~Ρ'k

 (k=1,2,..,Nz)とすれば,(3-4)のPoisson方程式は,

  

 k+1-2kk-1+Ωkk(k=2,3,..,Nz-1) (3-18a),

 21+Ω11 (3-18b),

 -NzNz-1+ΩNzNz (3-18c),

 

 すなわち,  

 Pi,j,k+1+{(hz/hy)2σj+(hz/hx)2λi-2}Pi,j,k+Pi,j,k-1

 =Ri,j,k   (k=2,3,.,Nz-1),

 

 Pi,j,2+{(hz/hy)2σj+(hz/hx)2λi -1}Pi,j,1

 =Ri,j,1,

 

  {(hz/hy)2σj+(hz/hx)2λi -1}Pi,j,Nz+Pi,j,Nz-1

 =Ri,j,Nz (i=1,2,.,Nx;j=1,2,.,Ny) 

 

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

  

 これはi,jを固定したとき,Nz次の三重対角行列の係数を持った

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

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

 解けます。

 

 1次元の連立1次方程式を解いてkが得られれば,

 それに対して,Φk=Q~kとすれば,最終的な解Φk,

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

 

 なお,Nx×Ny次の直交行列Q~=S~NyT~Nyの陽な形を書き下すと

 (Q~)ij,mn=Σk=1NxΣl=1Ny(S~Ny)ij,kl(T~Ny)kl,mn

 l(TNx)i,m(TNy)j,n となります。

 

 ただし,m(T)lm=(1/N)1/2 (m=1),

 (T)lm=(2/N)1/2cos{(2l-1)(m-1)π/(2N)}

 (m=2,3,..,N)です。...(3-20)

 

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

 

 第1行目が(-1,1,0,..,0),第k行目(k=2,3,..,N-1)が

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

 N次の三重対角行列をAとしすると,

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

 N個の異なる固有値λ=λm=-4sin2{(m-1)π/(2N)}

  (m=1,2,..,N)が得られます。

 

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

 θ=θm=(m-1)π/N (λ=λm=-4sin2m/2) )として,

 t(x1,x2,..,xN)の成分xkが,

 xk=a[exp(ikθ)+exp{-i(k-1)θ}}

 =2aexp(-iθ/2)cos{(2k-1)θ/2}なる形で与えられること

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

 

 そこで,λmに属する固有ベクトルを改めてmt(x1m,x2m,..,xNm)

 とおけば,xlm=cmcos{(2l-1)(m-1)π/(2N)}です。

 

 そして,固有値λ12,..,λN はすべて異なるので,

 m≠nなら実内積(mn)は直交します,つまり(mn)=tmn=0

 を満たします。

 

 一方,(mm)=tmm=cm2Σl=1Ncos2{(2l-1)(m-1)π/(2N)}

 =Ncm2/2 (m=2,3,..,Nのとき),Nc12 (m=1のとき)ですから,

 

 xlm=(2/N)1/2cos{(2l-1)(m-1)π/(2N)}(m=2,3,..,N),

 (1/N)1/2 (m=1)とおけば,(n)=tmn=δmnと正規直交化

 されます。

 

 そこで,N×N行列TNで成分が(TN)lm=xlmとなるものを与えると,

 このTN は1つの直交行列で,Aは直交変換tNATN=Λ(Λは対角

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

 

 (以上)

 

参考文献: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ランキングへ ← クリックして投票してください。(1クリック=1投票です。1人1日1投票しかできません。)

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (0) | トラックバック (0)

2007年5月19日 (土)

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

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

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

 

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

 

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

 

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

 

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

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

(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.

 

1.序論

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

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

 

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

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

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

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

 

   x(x1,x2,...,xn)を未知数とし,三重対角行列を係数とした

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

 

 a11+c12=y1,

 bkk-1+akk+ckk+1=yk (k=2,3,..,n-1),

 bnn-1+ann=yn

 

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

 

これを行列で表現すると,定数項を(y1,y2,..,yn)とし,

 

n次正方行列Aを,第1行目が(a1,c1,0,...,0),第k行目

(k=2,3,...,n-1)が(0,..,bk,ak,ck,0,..,0),第n行目

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

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

=Aと書けます。

 

  

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

 

1=a1,e1=y1として,dk+1=ak+1-ckk+1/dk,,

k+1=yk+1-ekk+1/dk(k=1,2,...,n-1)と順次計算して,

 

その後xn=en/dn,xn-k=(en-k-cn-kn-k+1)/dn-kとする,

 

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

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

A=QR(Rは上三角行列),Qは直交行列)とAをQR分解して,

=Q-1(e1,e2,...,en)とします。

 

Rをk行目(k=1,2,...,n-1)が(0,..,dk,ck,0,..,0),

第n行目が(0,...,0,dn)(ただしd1=1とする)の上三角行列

とすると,

 

A=QR,かつ=Qより,Qは第1行目が(1,0,...,0),

第k行目(k=2,3,..,n)が(0,..,bk/dk-1,1,0,..,0)で

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

そして1次元のPoisson方程式は,φ(x)を未知関数,ρ(x)を与えられた周知の関数とすると,d2φ/dx2=ρ(x)で与えられます。

 

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

k+1-2φk+φk-1)/h2=ρk (k=1,2,...,n)となります。

 

ただし,k=1のときはφ0を,k=nのときはφn+1を境界条件によって

与える必要があります。

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

1=(φ1-φ0)/h,f2=(φn+1-φn)/hを与える境界条件

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

φn+1を決めます。

 

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

[(φ2-φ1)/h-f1]/h=ρ1より,(φ2-φ1)/h2=ρ1+f1/h

が得られます。

 

そして,(-φn+φn-1)/h2=ρn-f2/hです。

  

ρ1をρ1+f1/hで,ρnをρn-f2/hで置き換えた列ベクトルを

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

このときn次の正方行列h2Aは第1行目が(-1,1,0,.. ,0),

第k行目(k=2,3,...,n-1)が(0,..,1,-2,1,0,..,0),

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

 

ところが,実はこのA,あるいはh2Aは固有値の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                 (注1終) ※

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

 矩形領域での2次元Poisson方程式は(∂2φ/∂x2)+(∂2φ/∂y2)=ρ(x,y) (a≦x≦b,c≦y≦d) ...(2-1)で与えられます。

 Neumann問題での境界条件はf(1),f(2),g(1),g(2)を既知関数として

 x=aで∂φ/∂x=f(1)(y),x=bで∂φ/∂x=f(2)(y),

 y=cで∂φ/∂y=g(1)(x),y=dで∂φ/∂y=g(2)(x)

  

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

 まず,2次元のメッシュの個数がNx×Nyであるとして差分hx,hy

 与えます:x0=a,xNx=b;xi-xi-1=(b-a)/Nx=hx

 (i=1,2,..,Nx),y0=c,yNy=d;yj-yj-1=(d-c)/Ny=hy

 (j=1,2,..,Ny) です。

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

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

 

 [(φi+1,j-2φi,j+φi-1,j)/hx2]

+[(φi,j+1-2φi,j+φi,j-1)/hy2]=ρi,j

 (i=1,2,...,Nx;j=1,2,...,Ny) ...(2-3)となります。

ただし,i=1のときにはφ0,j,j=1のときにはφi,0,i=Nxのとき

にはφNx+1,j,j=Nのときにはφi,Ny+1の境界値が与えられなけれ

ばなりません。

 

これらは境界条件(2-2)で,f(1)(y)をf(1)j,f(2)(y)をf(2)j,g(1)(x)をg(1)i,g(2)(x)をg(2)iと離散化し,以下のように書き直すことで与えられます:

 

φ0,j=φ1,j-hx(1)jNx+1,j=φNx,j+h(2)j,

φi,0=φi,1-hy(1)i,

φi,Ny+1=φi,Ny+hy(2)i

(i=1,2,..,Nx;j=1,2,..,Ny)...(2-4)

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

まとめると,

  

φj+1-2φjφj-1+(hy/hx)2φjρ'j

(j=2,3,..,Ny-1) ...(2-5a),

φ2φ1+(hy/hx)2φ1ρ'...(2-5b),

φNyφNy-1+(hy/hx)2φNy=ρ'Ny ...(2-5c)

 

となります。

ここに,φjt1,j2,j,..,φNx,j) (j=1,2,..,Ny),

ρ'jt1,j+f(1)j/hx2,j,..,ρNx,j-f(2)j/hx)

(j=2,3,..,Ny-1),

 

ρ'1t1,1+f(1)1/hx+g(1)1/hy ,ρ2,1+g(1)2/hy,..,ρNx,j

-f(2)j/hx+g(1)Nx/hy),

ρ'Nyt1,Ny+f(1)Ny /hx-g(2)1/hy2,Ny-g(2)2/hy..,ρNx,Ny-f(2)Ny/hx-g(2)Nx/hy)   ...(2-6)です。

また,Aは第1行目が(-1,1,0,..,0),第k行目(k=2,3,..,Nx-1)

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

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

Aの固有値方程式|λE-A|=0(Eは単位行列)は解析的に解くことができて,その固有値λ12,..,λNxは,

λm=-4sin2{(m-1)π/(2Nx)}  ...(2-8) で与えられます。

  (注2):以下,三重対角行列Aの固有値方程式|λE-A|=0 を具体的に解いてみます。

  

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

 

固有値方程式|λE-A|=0 はλについてn次の代数方程式なので,正攻法で解くのはかなり困難です。

 

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

 

k+1-(λ+2)xk+xk-1=0 (k=2,3,..,n-1)ですが,

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

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

ξ2-(λ+2)ξ+1=0 が得られます。

 

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

2cosθ=λ+2,すなわち固有値λはλ=2(cosθ-1)=-4sin2(θ/2)

と書くことができます。

 

そして,ξ=exp(iθ)がξ+1/ξ=λ+2の解なら,ξ-1=exp(-iθ)もそうですから,xk+1-(λ+2)xk+xk-1=0 の一般解はa,bを任意定数としてxk=aexp(ikθ)+bexp(-ikθ)と表わすことができます。

 

ここで,この一般解が境界条件;x2-x1=λx1,

λ=exp(iθ)+exp(-iθ)-2 に適合する条件を求めると, 

b=aexp(-iθ)となりますから,特にa=1,b=exp(iθ)

として,xk=exp(ikθ)+exp{-i(k-1)θ}とします。

 

これを-xn+xn-1=λxnに代入すると,sin(nθ)=0 を得ます。

 

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

θ=θ=(m-1)π/n (m=1,2,..,n)とすることができ

ます。

 

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

λ=λ=2(cosθ-1)=-4sin2/2)が,

固有値方程式|λE-A|=0 のn個の根を示していると

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

 Aは実対称行列で固有値はすべて異なる値なので,直交行列Tx

 (すなわち,tNxNx=TNxtNx=E)が存在して,

 

 AはΛNxtNxATNxNxは対角成分がλ12,..,λNx

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

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

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

 Tx(2/Nx)1/2(1,2,3..,Nx);

 

 1t(1/√2,1/√2,1/√2,..,1/√2),

 2=(cos{π/(2Nx)},cos{3π/(2Nx)},cos{5π/(2Nx)},..

 ,cos{(2Nx-1)π/(2Nx)}),

 3=(cos{2π/(2Nx)},cos{6π/(2Nx)},cos{10π/(2Nx)},..,

 cos{2(2Nx-1)π/(2Nx)}),..,

 

 Nx=(cos{(Nx-1)π/(2Nx)},cos{3(Nx-1)π/(2Nx)},

 cos{5(Nx-1)π/(2Nx)},..,cos{(2Nx-1)(Nx-1)π/(2Nx)})

 ...(2-10a),

 

 すなわち,

 (TNx)lm=(1/Nx)1/2 (m=1のとき),

 (TNx)lm=(2/Nx)1/2cos{(2l-1)(m-1)π/(2Nx)} 

 (m=2,3,..,Nxのとき) ...(2-10b)

 です。

 そこで,jtxφj,jtNxρ'j (j=1,2,..,Ny) ...(2-11)

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

 

j+1-2jj-1+(hy/hx)2ΛNxjj (j=2,3,..,Ny-1)

 ...(2-12a),

 21+(hy/hx)2ΛNx11 ...(2-12b),

 -NyNy-1+(hy/hx)2ΛNxNy=Ny ...(2-12c),

 

 すなわち,Pi,j+1+{(hy/hx)2λ-2}Pi,j+Pi,j-1=Ri,j

 (j=2,3,..,Ny-1),

 

 Pi,2+{(hy/hx)2λ-1}Pi,1=Ri,1,

 {(hy/hx)2λ-1}Pi,Ny+Pi,Ny-1=Ri,1Ny ...(2-13)

 

 となります。

 これらはi=1,2,...,Nxの各々を固定すると,j=1,2,...,Ny

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

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

 

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

を,'it(Pi,1,Pi,2,..,Pi,Ny),右辺の既知関数縦ベクトル

'it(Ri,1,Ri,2,..,Ri,Ny)とおき,行列A'i

 

第1行目が({(hy/hx)2λ-1},1,0,..,0),

第k行目(k=2,3,..,Nx-1)が,

(0,..,1,{(hy/hx)2λ-2},1,0,..,0),

第Nx行目が(0,,.,0,1,{(hy/hx)2λ-1})

 

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

未知数が'iの行列方程式; A'i'i'i (i=1,2,..,Nx) 

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

 

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

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

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

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

 

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

 以上の過程で求めたjに対してφj≡Txjによってφjを求めれば

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

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

  

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

 へと続きます。

  

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

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

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

 

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

 

 参考文献: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) 他

| | コメント (0) | トラックバック (0)

2006年11月11日 (土)

折り返しノイズ

前記事の続きです。

 

前記事ではf(t)のFourier変換:F(ω)≡(f)において,|ω|>ω1のときはF(ω)=0 であるという理想的なケースを考えました。

 

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

 

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

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

理想的なケースなら,(fT)=(2π/T){-∞(ω-2nπ/T)}の1つ1つのF(ω)の信号(ω-2nπ/T)でのゼロでない帯域,-ω1<ω-2nπ/T<ω1は決して重なることはないです。

 

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

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

 

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

 

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

 

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

 

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

 

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

  

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

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (5) | トラックバック (0)

2006年11月10日 (金)

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

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

 

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

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

 

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

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

(f)≡F(ω)≡∫-∞f(x)exp(-iωx)dxで定義します。

 

 このとき,逆変換は,f(x)=[1/(2π)]-∞(ω) exp(iωx)dωで与えられます。

また,周期がTの周期的デルタ関数:δT(x)を,

δT(x)≡∑-∞δ(x-nT)で定義します。

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

 

すなわち,δT(x)は周期Tの周期関数なので,これをFourier級数に展開すると,δT(x)=(1/T)∑-∞exp[i(2nπx/T)]となります。

 

そこで,(δT)=(1/T)-∞-∞[exp{-ix(ω-2nπ/T)}]dx(2π/T)-∞δ(ω-2nπ/T)と書けるからですね。

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

すなわち,(ω)≡(f)において,|ω|>ω1のときはF(ω)=0

であるとします。

 

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

 

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

ここで,ある時間的周期Tごとにこの信号(t)をサンプリングしたと考えると,サンプリングした信号をf(nT)として,全てのサンプル信号はfT(t)≡-∞(nT)δ(t-nT)=f(t)δT(t)と表現できます。

これをFourier変換すると,(fT(t))=[fδT]

[f(t)]*[δT(t)]となります。

 

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

 

そこで,(fT)=(2π/T)F(ω)*{-∞δ(ω-2nπ/T)}

(2π/T){-∞(ω-2nπ/T)}となります。

もしも,T(t)のFourier変換(fT)の周期(2π/T)が,周波数の幅である2ω1よりも大きい:(2π/T)>2ω1 or T<(π/ω1)なら,

 

1つ1つの帯域幅を持つ信号:(ω-2nπ/T)のゼロでない帯域:-ω1<(ω-2nπ/T)<ω1は決して重なることがないわけです。

  

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

 

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

具体的には,(fT)は周期(2π/T)の周期関数なのでFourier級数に展開できて,(fT)=(2π/T){-∞(ω-2nπ/T)}=-∞n exp(inTω);cn =∫-π/Tπ/TF(ω)exp(-inTω)dω=Tf(-nT)です。

 

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

そこで,(fT)=-∞(nT)exp(-inTω)であり,左辺は -π/T≦ω≦π/TではF(ω)に等しいので,

 

結局のところ,帯域内ではF(ω)=∑-∞(nT)exp(-inTω)と展開されます。

 

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

f(t)=∑-∞[f(nT)sin{π(t/T-n)}/{π(t/T-n)}]

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

 

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

 

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

 

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

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

 

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

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

 

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

 

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

 

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

 

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

 

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

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

にほんブログ村 科学ブログへクリックして投票してください。(ブログ村科学ブログランキング)

にほんブログ村 トラコミュ 物理学へ
    物理学

| | コメント (0) | トラックバック (0)

2006年9月30日 (土)

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

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

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

 

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

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

 

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

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

 

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

 それ故,1秒当たりの空気の流出入量をk(3/s)とするとk=nV/3600 ですね。

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

 

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

流入量はkC0Δtですから,部屋内のその量の時刻変化は,

V(C(t+Δt)-C(t))=k[C0-C(t)]Δtなる式で

与えられます。

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

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

dC/dt=a(C0-C)+Q/V と表わすことができます。

 

 ここで,a≡k/Vです。

 

 B=aC0+Q/VとおけばdC/dt=-aC+Bとなりますが,

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

 C(t)=(C( 0 )-a-1B)exp(-at)+a-1B 

 なる解が得られます。

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

 これを一般化して各部屋の容積がVi(i=1,2,..,m)のm個の部屋を有する住宅があって,これらの部屋の平均濃度が(m+1)次元の列ベクトルで(t)=t(C0,C1,C2,..,Cm)で表わされるとします。

 

 ここで,Ckはk番目の部屋における物質濃度です。

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

そして各部屋には排出強度Qiの排出源があり"部屋iと部屋jの間の換気量=先の方程式のkに相当する値"をkij=kji(i≠j;i,j=0,1,2,..m)とします。

  

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

dCi/dt=∑j=0mij(Cj-Ci)+Qi/Vi(i=0,1,2,..,m)

となります。ただしaij≡(kij/Vi)です。

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

 t(Q0/V0,Q1/V12/V2,..,Qm/Vm)で定義します。

 ただしQ0/V0=0 としておきます。

 

 また,λi≡∑j=0mijと定義し,行列Aを成分;(aij-λiδij)によって定義すれば,濃度を求めるボックスモデルの微分方程式:(m+1)変数の定数係数連立線形非同次微分方程式は,

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

 d/dt=-AB と表わされます。  

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

  

 初期値(0)に対し,(t)={(0)-A-1}exp(-At)+A-1

 と書くことができます。

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

 従来の慣用的な方法であれば,Aの最大(m+1)個の固有値:λi(i=0,1,2,..m)を全て求め,(t)の各成分を具体的にexp(-λit)の

1次結合で表わします。

 

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

 形式解:(t)={(0)-A-1}exp(-At)+A-1で,

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

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

 

 そこで,左辺のexp(-At)を右辺の∑ν=0のν=0,1,2,..の最初の数項で近似評価することによって力技で数値的に(t)の近似解を求めることができるわけです

 

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

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

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

| | コメント (0) | トラックバック (0)

2006年9月22日 (金)

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

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

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

 

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

 

 このときの拡散方程式の解である汚染物質濃度を,垂直距離x(m)と高さz(m)の関数と考えてC(x,z)で表わすと,これはAをQに比例するあるパラメータとして,C(x,z)=Ax-sexp(-Bzp/x)なる形の式(Robertsの式)で表わされる,ことがわかります。

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

   

 

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

 

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

そのパイプ状の発生源をy軸としたとき,n個の観測点座標:(xi,zi)(i=1,2,..n)において,風向が直角に近い環境のときの実測濃度:ciと先に設定した低煙源拡散式による計算値:Ci=C(xi,zi)とを比較することによって,逆にその式のパラメータA,s,B,pを推定します。

 具体的には,C(x,z)=Ax-sexp(-Bzp/x)の右辺をf(A,s,B,p;x,z)と書いて,i=f(A,s,B,p;xi,zi)として,この計算値と実測値ciの誤差の二乗和:S(A,s,B,p)≡∑(Ci-ci)2を最小にするパラメータA,s,B,pを求めるわけです。

 

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

 ここでは,より一般的に未知パラメータがr個あるとし,それらを順にA1, A2,..,Arとします。

 

 上のケースではr=4で,A1=A,A2=s,A3=B,A4=pです。

 

 そして,Ci=C(xi,zi)=f(A1,A2,..,Ar;xi,zi)とし,S(A1,A2,..,Ar)≡∑(Ci-ci)2と定義するわけです。

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

 

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

 i(A1,A2,..,Ar)≡∂S/∂Ai=2∑(Ck-ck)∂Ck/∂Ai=0 (i=1,2,..,r) で与えられます。

 

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

 

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

 具体的には,i(A1,A2,..Ar)=0 (i=1,2,..r)を線形近似することにより多変数ニュートン法(Newtonian method)を実行します。

 

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

 

 すなわち,初期値としてAj=Aj0 適当に与えた後に,

 

 0 = i(A1,A2,..Ar)

  =i(A10,A20,..,Ar0)+∑(∂Fi/∂Aj)|A=A0(j-Aj0)

 

 と近似します。

これを行列形式で書くため,Dという行列をD=(dij)≡(∂Fi/∂Aj)で定義し,特にD0(dij0)≡(∂Fi/∂Aj)|A=A0とします。

 

一方,列ベクトルA≡t(A1,A2,..,Ar)で定義し,特に

0t(A10,A20,..,Ar0)とします。

 

さらに,Fi0i(A10,A20,..,Ar0)を成分とする同様な列ベクトルを0と定義します。0t(F10,F20,..,Fr0)です。

 

こうすれば,先の線形近似は 00+D0(0)と表わされます。

 

これを単純に解くと,00-10と近似されることになります。

ここで0-10の逆行列です。

得られた00-10を,改めて0として代入してD0-10を計算し,00-10収束するまで繰り返します。

 

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

 

具体的には,100-10,211-11,..,m+1mm-1mを用います。ただし,Dk(dijk)≡(∂Fi/∂Aj)|A=Akです。

 

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

 

こうすれば,m→ ∞ でのmt(A1m,A2m,..,Arm)の極限値としてi(A1,A2,..,Ar)=0 (i=1,2,...r)を満たすt(A1,A2,..,Ar)が得られるはずです。

さらに,実際には収束を速くするために加速係数:wを与えて,m+1m-wm-1mとする方法などを用いた方がいいと思います。

  

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

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

 

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

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

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

| | コメント (2) | トラックバック (0)

より以前の記事一覧

その他のカテゴリー

001. 目次 | 002. 募金・ボランティア | 003. 日記・回想 | 004 訃報 | 005. 心身・思想・哲学 | 006. 社会・経済・政治 | 007. 病気(診察・薬) | 008. 恋愛・異性 | 009 宗教・神話 | 010 歴史(日本,世界) | 011. 将棋 | 012. TV(ニュース・ドラマ) | 013 スポーツ(ニュース・イベント) | 014 ノン・フィクション | 015 小説・詩・評論 | 016 漫画・劇画・アニメ | 017 演劇・映画・舞踊 | 018 音楽(日本・西洋・他) | 019 タレント(俳優・お笑い) | 020 ミュージシャン | 021 アイドル・ヒーロー | 022 創作 | 023 シャレ・ギャグ等 | 024 競馬・toto・賭け事 | 025 ファッション・風俗 | 100. 物理学一般 | 101 教育・学校(物理) | 102. 力学・解析力学 | 103. 電磁気学・光学 | 104. 熱力学・統計力学 | 105. 相対性理論 | 106. 星・ブラックホール・一般相対性 | 107. 重力・宇宙・一般相対性 | 108. 連続体・流体力学 | 109. 物性物理 | 110. 複雑系・確率過程・非線型・非平衡 | 111. 量子論 | 112. 原子・分子物理 | 113. 原子核物理 | 114 . 場理論・QED | 115. 素粒子論 | 116. 弦理論 | 118. 観測問題・量子もつれ | 119. 電気回路 | 200. 問題・解答 | 201. 自然科学一般 | 202. 気象・地学・環境 | 203. 生物学・生理学・生化学 | 204. 経済学(ミクロ・マクロ・マルクス) | 300 数学一般・算数 | 301. 集合・位相 | 302. 論理学・数学基礎論 | 303. 代数学・数論 | 304. 解析学 | 305. 複素数・複素関数論 | 306. 線型代数学 | 307. 幾何学(トポロジー・他) | 308. 微分方程式 | 309. 確率・統計 | 310. 関数解析・超関数 | 311 .数値計算・調和解析・離散数学 | 312. 公式・特殊関数 | 501. 商用宣伝・アフィリエイト