« 2007年6月 | トップページ | 2007年8月 »

2007年7月

2007年7月31日 (火)

北朝鮮で選挙

 北朝鮮で選挙があって投票率が99.82%だというので驚いた。100%じゃないのかよ。。。こりゃ大変だ。投票しなかった0.18%は刑務所行きか,収容所行きかと心配ですね。

 ヒトラーがいたナチス時代のドイツも選挙の投票率は100%だったと思ったが違ったかな。

 まあ,病気とか天災とか不可抗力と認められたら許されるのかな?日本でよかったですねえ。丸川さん。

 日本もマスコミを中心に選挙に行け,投票しろと結構うるさいけど,投票は権利であって義務じゃあない。まあ,時の政権にとっては投票率が低いほうが有利で野党にとっては逆らしい。

 そして大方の世論操作では「投票しない人間に政府や公共団体や代議員の政治に意見を言う資格はない。」などのキャンペーンがなされているらしいが,そんなことはない,投票しようがしまいが文句を言う資格はある。

 それに何年に数回かの選挙やその投票だけが政治活動であって,しかも何万票分の一,何十万票分の一じゃ,いかにも淋しい。

 かといって自分で議員というボランティアを一生の仕事としようというのも大変ですしね。

 しかも日本の現実じゃボランティアというには程遠いから,これから税金に巣くっている官僚などを本来の公僕にするという官僚のボランティア化も含めて体制変革から手がけるといって政党のようなものを作って政権を取ってから,シコシコとやろうとしても大変です。

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

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

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

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

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

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

遅い粘性流(5)(Stokes近似)

遅い粘性流の項目を終わりにするに当たって1982年5月の覚書きノートから「I.Imai(今井功)の方法(1972)によるストークス(Stokes)方程式の一般解」を紹介します。

まず,非圧縮性流体の定常流のストークス方程式と連続方程式はη△=gradp (1),div=0 (2)と表わされます。

 

そこで,今 △φ0 を満足するあるベクトル場:φ (3)によってp=2ηdivφ (4)と表わされたとします。

 

このとき(1),(2)の1つの解1は次式で与えられます。すなわち,1grad(rφ)-2φ (5) です。これが実際に解であることは1を直接(1),(2)に代入することで確かめられます。

次に一般解をとして,'≡1とおけば,'=0 ,div'=0 となります。第2の式によって,'=rotなるが存在します。

 

そして第1の式から,△'=rot(△)=0 なので,ある関数χが存在して△=gradχと書けます。

 

ここで,△Ψ=-χなるΨによって,ゲージ変換:'=+gradΨを施せば △'=0 とできることがわかります。

以上から,一般解1rot(6),△=0 (7)と表わすことができることがわかります。

 

そこで,χ≡-/2, △Φ=0 (8)なるχ,Φを用いてφ'=φ+gradΦ+rotχ(9)とおけば,△φ'=△φ=0 ,divφ'=divφ(10)が成立し,grad(rφ')-2φ'=1grad{gradΦ-2Φ+rot(-/2)}+rot(11)と書けます。

したがって,gradΦ-2Φ-(/2)rot=一定となるような調和関数Φが存在すれば,φ'1grad(rφ)-2φ(15)のφの代わりに用いて=grad(rφ')-2φ'とすることで,一般解が与えらます。

今,Φ=ΣlΦl,=Σll (12)と,Φ,をl次の体球関数Φl,l(調和関数であってx,y,zのl次の同次式):Φl=Σmlmllm(θ,φ),l=Σmlmllm(θ,φ))で展開すれば,(l-2)Φl(/2)rotlconst.δl0(13)が得られます。

そこで,20 でもrot20 のケースには,l≠2なるlに対するΦlとしてΦlrotl/{2(l-2)}と取れば,求めるΦ=ΣlΦlが存在することがわかります。

 

すなわち,Φ=(1/2)Σl≠2[rotl/(l-2)] (14) です。

さらにrot20 のケースでも(9)とは別のφ'の選び方で=grad(rφ')-2φ' (5)'の形に帰着できることが予想されます。

今井功 著「流体力学」(1973)(裳華房)

 

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

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

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

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

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

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

2007年7月30日 (月)

小田実逝く

 38年前,私が19歳のときに社会運動を始めた機会となり,最初に影響を受けた元ベ平連=「ベトナムに平和を!市民連合」代表の巨星=小田実(おだまこと)氏が75歳で亡くなったらしいです。

   合掌!

 まあ,市民連合の市民というのは"ブルジョア"ですが,そういう意味の「市民連合」ではなかったと思います。

 同じ頃,市民運動を名乗って東京工大でビラを配っていたらしい「菅直人」という俗物とは対極に位置していますね。

      。

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

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

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

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

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

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

遅い粘性流(4)(Oseen近似)

1982年6月の覚書きノートから抜粋して,1928年のゴールドスタイン(S.Goldstein)の方法により,球を過ぎる一様流に対する定常オセーン近似の解について記述します。

半径aの小球のまわりの線型粘性流の速度ベクトル()で表わし,lim|r|→∞=Uとします。

 

=Uとおけば,に対するオセーン近似(Oseen近似)の方程式と連続の方程式は,U(∂/∂x)=-gradp/ρ+η△/ρ (1),div=0 (2)と表わされます。ρは密度,ηは第一粘性率です。

小球の中心を原点とし,r=||とおけばに対する境界条件は,=-Uat r=a (3),0 at r=∞ (4),p=p(一定) at r=∞ (5)と書けます。

(1),(2)より,△p=0 (6)なのでp=ρU(∂Φ/∂x)(あるいは,p=ρU(∂Φ/∂x)+p(7),かつ△Φ=0 (8)を満たす関数Φ()が存在します。

なぜなら,一般にp=ρU(∂Φ/∂x)なら,∂(△Φ)/∂x=0 より,Fをy,zの任意関数として△Φ=F(y,z)と書けます。

 

特にΦ≡{1/(ρU)}∫-∞x(x',y,z)dx'と置くと,r→ ∞ に対してp → p(一定)より∂p/∂x → 0 ですから,ρU△Φ=(∂p/∂x)+∫-∞x [(∂2/∂y2)+(∂2/∂z2)]=(∂p/∂x)-[(∂p(x',y,z)/∂x')]-∞x=0 が成立するからです。

ここで,特に≡-gradΦ (9)とおけば(7),(8),(9)は明らかに(1),(2)の1つの解であることがわかります。

 

そこで,一般解を求めるために'を≡-gradΦ+'(10)で定義します。また,k=ρU/(2η) (11)とおきます。

このとき,(1),(2)は,それぞれ[△-2k(∂/∂x)]'=0 (12),divu'=0 (13)となります。

 

渦度ωω≡rot=rotu' (14)で定義すると,[△-2k(∂/∂x)]ω=0 (15)が成立します。

ところで,x軸(θ=0)の周りで流れは完全に対称なので,この軸を含む任意の平面に関して流速ベクトルは左右対称であり,その平面を横切るような流速成分は存在しません。

 

したがって,uφ0 であり渦線は極軸(θ=0)のまわりの円環を形成することがわかります。

 

すなわち,ωは上述の平面に垂直な成分のみを持つのでω=0 ,つまりω0 (16) です。

故に,div(rot)=0 より(∂ωy /∂y)+(∂ωz /∂z)=0 (17)を得ます。そこで,xと(y2+z2)1/2のみ,つまりrcosθとrsinθのみの関数χ(r,θ)が存在してωy =-∂χ /∂z,ωz =-∂χ /∂y (18)と書けます。すなわち,ω=-rot(χ)です。

(15)より,∂[{△-2k(∂/∂x)}χ]/∂z=0,∂[{△-2k(∂/∂x)}χ]/∂y=0 なので,{△-2k(∂/∂x)}χ=g(x)と書けます。

 

ここで,df/dx=-e2kx∫g(x)e-2kxdx;g(x)≠0 なるf(x)に対してχ'≡χ+f(x)と置けば{△-2k(∂/∂x)}χ'=0 とできるので,一般性を失うことなく{△-2k(∂/∂x)}χ=0 (19)としてよいことになります。

ところで,rotω=rot(rot')=grad(div')-△'=-△'ですが,一方,ω=-rot(χ)よりrotω=-grad{div(χ)}+(△χ)です。

 

これと,{△-2k(∂/∂x)}χ=0 (19)よりrotω=-∂(gradχ-2kχ)/∂x (20)を得ます。また,[△-2k(∂/∂x)]'=0 (12)よりrotω=-2k(∂'/∂x) (21)も成立します。

そこで,(20),(21)より∂{'-gradχ/(2k)+χ}/∂x=0 (22)なので,'≡ gradχ/(2k)-χ(y,z)と書けます。

 

ところで,div'=0,rot'=-rot(χ),[△-2k(∂/∂x)]'=0,{△-2k(∂/∂x)}χ=0 なので,div=0,rot=0,△=0 となります。

はy,zだけの関数なので,y,zだけの関数ψ(y,z)が存在して=-gradψ,△ψ=(∂2ψ/∂y2)+(∂2ψ/∂z2)=0 と書けます。

 

そして,u=-gradΦ+' (10)よりu=-grad(Φ+ψ)+gradχ/(2k)-χとして良いことになります。∂ψ/∂x=0 ですからp=ρU(∂Φ/∂x) (7)はp=ρU{∂(Φ+ψ)/∂x} (7)'と書けます。

また.△Φ=0 (8)も△(Φ+ψ)=0 (8)'と書けます。それ故,ベクトル=-gradψは解にとって本質的な量ではないので落としてもかまいません。

したがって,(22)より'= gradχ/(2k)-χ(23),あるいは=-gradΦ+gradχ/(2k)-χ(24)と書いてよいことになります。

 

ところで,χ'≡e-k xχ ⇔ χ≡ekxχ'(25)によってχ'を定義すれば{△-2k(∂/∂x)}χ=0 (19)によってχ'に対してヘルムホルツの方程式(△-k2)χ'=0 (26)が成立します。

次に解を極座標系(r,θ,φ)で表わすことを考えます。r,θ,φ方向単位ベクトルをそれぞれr,θ,φと置きます。

 

任意関数fについてdf=(gradf)d=(gradf)(drr+rdθθ+rsinθdφφ)なので,gradf=(∂f/∂r)r(1/r)(∂f/∂θ)θ{1/(rsinθ)}(∂f/∂φ)φです。 

また,i=e1(r,θ,φ)t(cosθ,-sinθ,0),j=e2(r,θ,φ)t(sinθcosφ,cosθcosφ,-sinφ),k=e3(r,θ,φ)t(sinθsinφ,cosθsinφ,cosφ)となります。

 

それ故,(i,j,)=(1,2,3)≡(r,θ,φ)Tで変換行列Tを定義すれば,t(ur,uθ,uφ)=Tt(u1,u2,u3)となります。

それ故,=-gradφ+gradχ/(2k)-χ(24)を極座標表示すると

 

r=-(∂Φ/∂r)+{1/(2k)}(∂χ/∂r)-χcosθ=-(∂Φ/∂r)+{ekx/(2k)}(∂χ'/∂r)-(ekx/2)χ'cosθ(27),

 

θ=-(1/r)(∂Φ/∂θ)+{1/(2kr)}(∂χ/∂θ)+χsinθ=-(1/r)(∂Φ/∂θ)+{ekx/(2kr)}(∂χ'/∂θ) +(ekx/2)χ'sinθ(28),uφ0 (29)となります。

また,Φ,χ'をr,θで表わすと次式のようになります。

 

Φ=Σn=0(An/rn+1)Pn(cosθ) (30),χ'=Σn=0[Bnχn(kr)Pn(cosθ)] (31)です。ここに,χn(z)≡(2n+1){π/(2z)}1/2n+1/2(z) (32)です。

境界条件(3),(4)は,次のようになります。ur=-Ucosθ,uθ=-Usinθ at r=a (33),ur=uθ0 at r=∞ (34)です。

(33)を用いると,半径aの球全体に働く抵抗をΦとρ,Uのみを使って表現することができます。

 

以下にそれを証明します。

ニュートン・ストークスの公式によれば応力テンソルPijはPij=-pδij2ηeij (35)で与えられます。

 

ここで,eijは歪み速度テンソルでeij(∂vi/∂xj+∂vj/∂xi)/2=(∂ui/∂xj+∂uj/∂xi)/2 (36)で定義されます。

 

そして,球面の外向き法線として,応力を()で表わすとはfi=Pijj (37)で与えられます。

球面の外向き法線方向は常にr方向に一致するので,(35),(37)よりfr=Prr=-p+2ηerr (38),fθ=Pθr=P2ηe (39),fφ=Pφr=P2ηe (40)となります。

直交曲線座標に関する定理を,線素dr,rdθ,rsinθdφに適用すると,err,e,eは次のように表わされます。

 

rr=∂ur/∂r (41),e[(1/r)(∂ur/∂θ)+r{∂(uθ/r)/∂r}]/2 (42),e[{1/(rsinθ)}(∂ur/∂φ)+r{∂(uφ/r)/∂r}]/2≡0 (43)です。よって,fφ=P0 です。

そして,抵抗は次式で与えられます。

 

r=a[Prrt(cosθ,sinθcosφ,sinθsinφ)+t(-sinθ,cosθcosφ,cosθsinφ)]dS (44) です。

 

r=Prrおよびfθ=Pは角変数φを含まないので,明らかに2=D30 です。

 

そして,ゼロではない唯一の抵抗成分D1を改めてDと書けば,D=2πa20πdθsinθ(Prrcosθ-Psinθ) (45)となります。

ここで,まず圧力pをφ,χで表現すると,p=ρU(∂Φ/∂x)=ρU{cosθ(∂Φ/∂r)-(sinθ/r)(∂Φ/∂θ)} (46)です。

また,(27),(41)よりerr=∂ur/∂r=-(∂2Φ/∂r2)+{1/(2k)}(∂2χ/∂r2)-(∂χ/∂r)cosθです。

 

△Φ=0 より,-(∂2Φ/∂r2)=(2/r)(∂Φ/∂r)+{1/(r2sinθ)}{sinθ(∂Φ/∂θ)}です。

 

また,{△-2k(∂/∂x)}χ=0 から,{1/(2k)}(∂2χ/∂r2)={-1/(kr)}(∂χ/∂r)-{1/(2kr2sinθ)}{sinθ(∂χ/∂θ)}+(∂χ/∂r)cosθ-(sinθ/r)(∂χ/∂θ)ですから,結局次式が得られます。

rr(-2/r)[-(∂Φ/∂r)+{1/(2k)}(∂χ/∂r)-χcosθ]-{1/(rsinθ)}(∂/∂θ)[(-1/r)(∂Φ/∂θ)+{1/(2kr)}(∂χ/∂θ)}+χsinθ],つまりerr=-2ur/r-{1/(rsinθ)}{∂(uθsinθ)/∂θ}=2Ucosθ/a-2Ucosθ/a≡0 (47)となります。

また,(27),(28),(42)よりe[(1/r)(∂ur/∂θ)+r{∂(uθ/r)/∂r}]/2={-uθ/r(∂uθ/∂r)+(1/r)(∂ur/∂θ)}/2は次のように変形されます。

ここで,ノートには詳細な計算過程が記されていますが,eについてもrr=∂ur/∂rの場合と同様なので,詳細は割愛します。

 

結局,e=ρU{(∂Φ/∂r)sinθ+(cosθ/r)(∂Φ/∂θ)} (48) となります。

(38),(39),および(47),(48)によりPrrcosθ-Psinθ=-ρU(∂Φ/∂r) (49)となります。

 

そこで,小球が流れの向きに受ける抵抗力の大きさは,D=2πa20πdθsinθ(Prrcosθ-Psinθ) (45)によりD=-2πa2-11(cosθ)[ρU(∂Φ/∂r)] (50)となります。

 

そして,Φ=Σn=0(An/rn+1)Pn(cosθ) (30)によれば∂Φ/∂r=-Σn=0{(n+1)Ann(cosθ)/rn+2}なので,Pn(cosθ)の直交性により,結局D=4πρUA0 (51)なる重要な式が得られました。

以下は,もうノートには記載されていなくて手元にはゴールドスタイン(S.Goldstein)の論文しかないので概略のみ記述します。

が境界条件:ur=-Ucosθ,uθ=-Usinθ at r=a (33),ur=uθ0 at r=∞ (34)を満たすという条件から次式を得ます。

 

Φ=Σn=0(An/rn+1)Pn(cosθ) (30),χ'≡e-k xχ (25),χ'=Σn=0[Bnχn(kr)Pn(cosθ)] (31) (ただしχn(z)≡(2n+1){π/(2z)}1/2n+1/2(z))における係数AnとBnの関係が得られ,A0{π/(4k2)}Σm=0{(2m+1)Bm}です。

 

つまり,k=ρU/(2η) (11),D=4πρUA0 (51)から,D=(2π2η/k)Σm=0{(2m+1)Bm},あるいは抵抗係数をCD≡D/(ρπa22/2)とおいて,CD{2π/(ka)2m=0{(2m+1)Bm /U}となります。

レイノルズ数RはR=ρUL/η=UL/ν=2ρUa/η=4kaとなるので,境界条件によってBnの満たすべき条件から,球Bessel関数χn(ka)を近似計算して抵抗係数CDを求め,(ka)あるいはRのべきで展開します。

 

それによって,この論文の結論となる近似式:CD(24/R)[1+3R/16-19R2/1280+71R3/20480-30179R4/34406400+122519R5/560742400+...] (52)が得られます。

しかし,先に指摘したようにオセーン近似は遅い流れに対する単なる近似であり,そして近似としても重大な欠点を持っています。

 

それ故,この展開式(52)はオセーン近似の範囲内では正しくてもナビエ・ストークス方程式に基づく厳密な抵抗係数とはRの高次の項で一致しません。

 

したがって,この式のRの高次の展開項はある意味では意味がないと考えられます。

ナビエ・ストークス方程式に基づく厳密解によれば,これと同一の流体系での抵抗係数CDとしてCD(24/R)[1+3R/16+(9R2/160)(logR+γ+2log2/3-323/360)+27R3logR/640+O(R3)] (53)なる表式が得られています。ここにγはEulerの定数γ=0.57721...です。

 

これを見ると,近似式(52)は第2項,つまりRの1次までしか厳密解を正しく表現していないことがわかります。

 

参考文献:Sydney.Goldstein,B.A.,Ph.D.St.John's College,Cambridge."Steady Flow of Viscous Fluid Spherical Obstacle at Small Reynolds Numbers"(Communicated by H.Jeffreys,F.R.S.-Received December31,1928) Proc.Roy.Soc.A 123,pp225-235,巽友正 著「流体力学」(1982)(培風館)

 

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

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

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

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

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

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

2007年7月29日 (日)

遅い粘性流(3)(Oseen近似)

非圧縮かつ密度が一様な単成分ニュートン流体の運動方程式はナビエ・ストークス方程式ρ(∂/∂t)+ρ(・grad)=-gradp+η△と"連続の方程式=質量保存の方程式"div=0 で与えられます。

 

ここでρは流体の密度,(,t)は流体の速度,pは圧力,ηは第一粘性率です。そして,特に∂/∂t≡0 で流速が時刻tによらないとし,無限遠において流れが一様=U(r=||→ ∞)となる流体系を考えます。

レイノルズ数(Reynolds number)Rが小さい遅い粘性流について,ストークス近似(Stokes近似)ではナビエ・ストークス方程式の右辺の粘性項と比較して"左辺第二項=慣性項"全体を無視しましたが,今回はより高次の近似を考えるため,一様流に対する摂動を'とします。

 

すなわち,=U'によって'を定義します。明らかに'=0 at ||→∞ です。

そこで'の2次以上を無視すると,定常ナビエ・ストークス方程式は U(∂'/∂x)=-gradp/ρ+η△'/ρ (1)となり,連続方程式はdiv'=0 (2) となります。

 

この近似をオセーン(Oseen)近似といいます。

以下,ストークス近似のときと同様,半径aの小球の中心を原点とし,そのまわりの遅い流れについての境界値問題を考察します。

に対する境界条件が0 at r=||=a;=U,p=p at r=∞ なので,v'に対する境界条件は'=-Uat r=a (3),'= 0 at r=∞ (4),p=p at r=∞ (5) となります。

しかしながら,(3)から明らかなように球の付近では,摂動の仮定であるU>>|'|という条件が全く成立しません。

 

したがって,そこではむしろストークス近似の方が妥当な近似になります。オセーン近似というのは実際にはこうした欠点があるということを承知の上での近似であることに留意しておくことは重要です。

以下では,ストークス近似と同様な方法で'の解を求めることを試みることにします。

 

まず,(1)の発散を取り,(2)を適用すると△p=0 (6) となります。したがってpは軸対称な調和関数です。x軸を極軸に取るとpは次式で与えられます。p=p+Σl=0(Al/rl+1)Pl(cosθ)(7) です。

一方,オセーン近似の方程式:(∂'/∂x)=-gradp/ρ+η△'/ρ (1)は次のように変換することによって非斉次のヘルムホルツ方程式に帰着します。

すなわち,'≡~ekx,k=ρU/(2η) (8)とおけば,∂'/∂x=ekx{(∂~/∂x)+k~},△v'=ekx{△~+2k(∂~/∂x)+k2~}となります。

 

それ故,△'-(ρU/η)(∂'/∂x)=ekx(△-k2)~ですから,結局,オセーン近似の方程式は,非斉次ヘルムホルツ方程式:(△-k2)~=(e-kx/η)gradp (9) になります。

 

一方,連続方程式:div=0 は~の式ではdiv~=-kv~x (10) と表現されます。

また, ~に対する境界条件は次のようになります。すなわち, ~=-Ue-kxat r=a (11),および~= 0 at r=∞ (12) です。

ヘルムホルツ演算子(△-k2)に対する第一種グリーン関数を次の2つの式を満足するG(,')によって定義します。

 

すなわち,(△-k2)G(,')=δ(') (13), G(,')=0 at r=a,∞ (14) です。また,対称性:G(',)=G(,')を要求します。

グリーンの公式によれば,∫(u△v-v△u)dV=∫[u(∂v/∂n)-v(∂u/∂n)]dSですから,∫[u(△-k2)v-v(△-k2)u]dV=∫[u(∂v/∂n)-v(∂u/∂n)]dSです。

 

故に,∫[~(△-k2)G-G(△-k2)~]dV=∫[~(∂G/∂n)-G(∂~/∂n)]dSが成立します。

グリーン関数の性質(13),(14)と対称性,そして~が(9)~(12)の解であることから次式が成立します。

 

~()=(1/η)∫r'=ar'=∞(,')e-kr'cosθ'grad'p(')d3'+U∫r'=∞-kr'cosθ'{∂G(,')/∂r'}dS'(15)です。

 次に,第一種グリーン関数G(,')を具体的に求めます。

よく知られているように,球対称であって,無限遠でゼロとなるようなヘルムホルツ演算子のグリーン関数は,湯川ポテンシャル-e-k|r-r'|/(4π|'|)で与えられます。

 

つまり,(△-k2)(e-k|r-r'|/|'|)=-4πδ(')です。

 

そこで(13),(14)を満足するグリーン関数は,この湯川ポテンシャルに斉次ヘルムホルツ方程式(△-k2)ψ()=0 の特殊解でr>aの領域で特異点を持たないものを加えることで得られるはずです。

そこで,ヘルムホルツ方程式:(△-k2)ψ()=0 を解きます。'を極軸に取り'に対する偏角をωとするとG(,')は極軸に関して対称であり湯川ポテンシャルもまたそうです。

 

そこで,ψ()も極軸に関して対称であるべきなのでψ()=ψ(r, ω)とおけば,△ψ=(1/r2)[(∂/∂r){r2(∂ψ/∂r)}]+{1/(r2sinω)}[(∂/∂ω){sinω(∂ψ/∂ω)}]となります。 

特に,ψ(r,ω)≡R(r)Ω(ω)とおいて(△-k2)ψ(r,ω)=0 の変数分離の解を求めます。

 

このとき,(△-k2)ψ=0 は(1/r2)[(d/dr){r2(dR/dr)}] Ω+(R/r2)(1/sinω)[(d/dω){sinω(dΩ/dω)}]-k2RΩ=0 となります。

 

それ故,(1/sinω)[(d/dω){sinω(dΩ/dω)}]+λΩ=0,(d/dr){r2(dR/dr)}-(k22+λ)R=0 が得られます。

この固有値λに対する固有値問題のω部分Ωについてのλに属する固有解の全てはΩλ(ω)=Ωl(ω)=Pl(cosω),λ=l(l+1) (ただしl=0,1,2...)で与えられます。

 

そこで,動径部分はd2l/dr2(2/r)(dRl/dr)-{k2+l(l+1)/r2}Rl0 を満たします。

これの解Rl(r)は線形常微分方程式の解を与える公式集によれば,Rl(r)=Jl+1/2(kr)/r1/2です。

 

右辺のJνはν次のベッセル関数であり,この表現は球ベッセル関数と呼ばれています。

 

あるいは,B,Cを任意定数としてRl(r)=BIl+1/2(kr)/r1/2+CKl+1/2(kr)/r1/2とも書けます。

 

ここにIν,Kνはそれぞれν次の第一種,第二種の変形ベッセル関数です。

これらのうちで,r→∞ で有限な解はKl+1/2(kr)/r1/2の定数倍で与えられるので,結局ψ(r,ω)=Σl=0{cll+1/2(kr)/r1/2}Pl(cosω)と表わすことができます。

 

かくして,G(,')=-e-k|-r'|/(4π|'|)+Σl=0{cll+1/2(kr)/r1/2}Pl(cosω)(r≧a)です。

ここで,以前のストークス近似でも記述したように=(r,θ,φ),'=(r',θ',φ')なら,cosω=cosθcosθ'+sinθsinθ'cos(φ-φ'),Pl(cosω)=Pl(cosθ)Pl(cosθ')+2Σm=1l{(l-m)!/(l+m)!}Pml(cosθ)Pml(cosθ')cos{m(φ-φ')}が成立します。

 

これらの公式とPl,Pmlの直交性を用いて,r=aでG(,')=0 となるような係数clを求めます。 

まず,e-k|-r'|/|'|を変形ベッセル関数とルジャンドル多項式の積の級数に展開すると,e-k|-r'|/|'|=(rr')-1/2Σl=0{(2l+1)Il+1/2(kr)Kl+1/2(kr')}Pl(cosω)(a≦r<r'),およびe-k|-r'|/|'|=(rr')-1/2Σl=0{(2l+1)Il+1/2(kr')Kl+1/2(kr)}Pl(cosω) (a≦r'<r) となります。

それ故,G(,')=[-1/{4π(rr')1/2}]Σl=0(2l+1)Pl(cosω){Il+1/2(kr)Kl+1/2(kr')-Il+1/2(ka)Kl+1/2(kr')Kl+1/2(kr)/Kl+1/2(ka)} (a≦r<r') (16-1),

 

および,G(,')=[-1/{4π(rr')1/2}]Σl=0(2l+1)Pl(cosω){Il+1/2(kr')Kl+1/2(kr)-Il+1/2(ka)Kl+1/2(kr')Kl+1/2(kr)/Kl+1/2(ka)} (a≦r'<r) (16-2) を得ます。

 

これは確かに,'について対称です。

求める第一種のグリーン関数が(16-1,2)式によって完全に求まったので,~()=(1/η)∫r'=ar'=∞(,')e-kr'cosθ'grad'p(')d3'+U∫r'=∞-kr'cosθ'{∂G(,')/∂r'}dS' (15) に従ってv~()を求めます。

まず,(15)の右辺第二項(剰余項)を計算します。

 

-kacosθ'{π/(2ka)}1/2Σl=0(-1)l(2l+1)Pl(cosθ')Il+1/2(ka)であり,一方[Il+1/2(kr)/r1/2}'=-(1/2)r3/2l+1/2(kr)+kr-1/2I'l+1/2(kr),[Kl+1/2(kr)/r1/2}'=-(1/2)r3/2l+1/2(kr)+kr-1/2'l+1/2(kr)です。

 

そして,I'ν(z)Kν(z)-Iν(z)K'ν(z)=1/zを用いると,{∂G(,')/∂r'}r'=a[-1/{4πa(ar)1/2}]Σl=0[(2l+1)Pl(cosω)Kl+1/2(kr)/Kl+1/2(ka)]なので,次式が得られます。

すなわち,U∫r'=∞-kr'cosθ'{∂G(,')/∂r'}dS'=a2U∫0dφ'∫-11(cosθ'){∂G(,')/∂r'} r'=aより,(15)の"右辺第二項(剰余項)"=U∫r'=∞-kr'cosθ'{∂G(,')/∂r'}dS'=-U{π/(2k)}1/2Σl=0[(-1)l(2l+1)Pl(cosθ){Il+1/2(ka)/Kl+1/2(ka)}{Kl+1/2(kr)/r1/2}] (17)を得ます。

ここで,計算においては三角関数とルジャンドル多項式の直交性を用いました。

(15)式の右辺第一項を求めるために,p=p+Σl=0(Al/rl+1)Pl(cosθ).(7)からgradpを求めます。

 

これは既に前記事で得られた表式と同じです。

 

∂p/∂x=-Σl=0{(l+1)Al/rl+2}Pl+1(cosθ)(18),∂p/∂y=-Σl=0(Al/rl+2)P1l+1(cosθ)cosφ(19),∂p/∂z=-Σl=0{(l+1)Al/rl+2}P1l+1(cosθ)sinφ (20) ですね。

~()を与える(15)式の右辺第一項の表式は(1/η)∫r'=ar'=∞(,')e-kr'cosθ'grad'p(')d3'=(1/η)∫0dφ'∫-11(cosθ')∫adr'r'2[{π/(2ka)}1/2Σl=0(-1)l(2l+1)Pl(cosθ')Il+1/2(ka)]G(,')grad'p(')です。

 

これに,(16)のG(,'),および(18)~(20)のgrad'p(')の表現式を代入します。

というところで当時は力尽きて未完となっています。

 

ここから先は当時会社の先輩で契約社員だったN大理工学部助手(当時)のM前(Mさき)さんに頼んで,N大の図書館からこの計算についての1928年のゴールドスタイン(Goldstein)の論文のコピ-を取ってきてもらって,フォローしています。

 

その内容についてのノートは1982年の6月という2年後に移行しています。そこで切りがいいので今日はここまでにします。

 

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

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

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

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

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

icon

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

2007年7月28日 (土)

遅い粘性流(2)(Stokes近似)

前記事の続きです。

次に(10),(11),(12)の∂p/∂x,∂p/∂y,∂p/∂zの表現を用いて()=∫r'=ar'=∞(,')gradp(')d3'+ηU∫r'=∞∂G(,')/∂r'dS' (8) の右辺第一項を求めます。

x成分は,vx-U(1-a/r)={1/(2η)}Σl=0{2(l+1)All+1(cosθ)/(2l+3)}[∫0r{r'l+1/rl+2-a2l+1/(rr')l+2}(1/r'l)dr'+∫r{rl+1/r'l+2-a2l+1/(rr')l+2}(1/r'l)dr']です。

 

この積分を実行すると,vx=U(1-a/r)+{1/(2η)}Σl=0{(l+1)All+1(cosθ)/(2l+1)}{(r2-a2)/rl+2} (15) となります。

同様に,vy{1/(2η)}Σl=0{(l+1)Al1l+1(cosθ)cosφ/(2l+1)}{(r2-a2)/rl+2}(16),およびz{1/(2η)}Σl=0{(l+1)Al1l+1(cosθ)sinφ/(2l+1)}{(r2-a2)/rl+2}(17)を得ます。

圧力pの一般形:p=p+Σl=0(Al/rl+1)Pl(cosθ) (9)と,今求めたばかりのに対する(15),(16),(17)の表式により,境界条件(4),(5)と△p=0 (3)を満たす圧力pとポアソン方程式:η△=gradp (1)の解として流速が得られたことになります。

しかし,これらが連続の方程式:div=∂vx/∂x+∂vy/∂y+∂vz/∂z=0 (2)を満たすかどうかはわかりません。

 

そこで,具体的に(15),(16),(17)のdivを計算し,それをゼロとおくことによって係数Alの満足すべき条件を求めます。

泥臭い計算によって,∂vx/∂x=aUcosθ/r2{1/(2η)}Σl=0(l+1)Al[{-lr2(l+2)a2}cosθPl+1(cosθ)+(r2-a2)sin2θ{dPl+1(cosθ)/d(cosθ)}]/{(2l+1)rl+3},

 

∂vy/∂y={1/(2η)}Σl=0l[{-lr2(l+2)a2}sinθcosφP1l+1(cosθ)+(r2-a2)(cosθsinθcosφ{dP1l+1(cosθ)/d(cosθ)}-(sinφ/sinθ)P1l+1(cosθ))]/{(2l+1)rl+3},

 

∂vz/∂z={1/(2η)}Σl=0l[{-lr2(l+2)a2}sinθsinφP1l+1(cosθ)+(r2-a2)(cosθsinθsinφ{dP1l+1(cosθ)/d(cosθ)}-(cosφ/sinθ)P1l+1(cosθ))]/{(2l+1)rl+3}

 

が得られます。

sinθ{dPl+1(cosθ)/d(cosθ)}=P1l+1(cosθ),sinθ{dP1l+1(cosθ)/d(cosθ)}=(cosθ/sinθ)P1l+1(cosθ)-(l+1)(l+2)Pl+1(cosθ)を用いて変形します。

 

すると,div=aUcosθ/r2(1/η)Σl=0l[sin2θ{dPl+1(cosθ)/d(cosθ)}+(l+1)cosθPl+1(cosθ)]/{(2l+1)rl+1}となります。

さらに,公式(1-z2){dPl(z)/dz}+lzPl(z)=lPl-1(z)によって,sin2θ{dPl+1(cosθ)/d(cosθ)}+(l+1)cosθPl+1(cosθ)=(l+1)Pl(cosθ)です。

 

結局,div=aUcosθ/r2(1/η)Σl=0(l+1)All(cosθ)/{(2l+1)rl+1}=0 (18)を得ます。

 

これが成立するためには,2Al/(3η)=-aUδl1,つまりl=-(3/2)ηaUδl1(19)となるのが必要十分です。

したがって,p=p-(3/2)ηaUcosθ/r2=p-(3/2)ηaUx/r3 (20)です。

 

また,u≡vx=U(1-a/r)-{aU/(4r)}(1-a2/r2)(3cos2θ-1)=U-{aU/(4r)}(3+a2/r2)-{3aUx2/(4r3)}(1-a2/r2)(21),v≡vy=-{3aU/(4r)}(1-a2/r2)cos2θsinθcosφ=-{3aUxy/(4r3)}(1-a2/r2) (22),w≡vz=-{3aU/(4r)}(1-a2/r2)cos2θsinθsinφ=-{3aUxz/(4r3)}(1-a2/r2) (23)も得られます。

かくして,方程式(1)~(5)は完全に解かれ,解は一意的であることも示されました。

方程式が解かれたので,この解を利用して流体が球に及ぼす力を求めることができます。

 

この力は球面上に働く応力(単位面積当たりに働く力;法線および接線応力)を小球面にわたって積分することで求めることができます。

まず,一般的な応力テンソルPijを考えます。運動がナビエ・ストークスの方程式に従う流体は,応力テンソルが次のニュートン・ストークスの法則に従う線形ニュートン流体です。

 

すなわち,応力テンソルPijはPij(-p+ζellij2η(eij-ellδij/3) (i,j=1,2,3はそれぞれx,y,z成分に対応) (24)で与えられます。

 

ここで,pは圧力,η,ζはそれぞれ第一粘性率,第二粘性率であり,eij(∂vi/∂xj+∂vj/∂xi)/2は歪み速度テンソルです。

非圧縮性流体では連続の方程式としてdiv=ell0 が成立するので,(24)はPij=-pδij+η(∂vi/∂xj+∂vj/∂xi) (25)となります。

 

応力と応力テンソルPijの関係はfi=Pijj(jはアインシュタインの規約に従うダミー添字)で与えられます。

 

また,は応力の働くべき面の外向き法線ベクトルで,小球面上では(nx,ny, nz)=/r=(x/r,y/r,z/r)r=a(cosθ,sinθcosφ,sinθsinφ)と表わされます。

小球全体に働く力は応力によって=∫dSで与えられます。ところが流れはx軸に関して対称(uはy,zの偶関数),v,wはy,zの奇関数)ですから,∫f2dS=∫f3dS=0 となるはずです。

 

実際,f2=P2jj,f3=P3jjなので,解からf2(-)=-f2(),f3(-)=-f3()あるいはP2j,3jjの偶関数であることが示されます。

まず,-pδ2j,-pδ3jですが,p=p-(3/2)ηaUx/r3は明らかにy,zの偶関数であり,よって-pδ2j,-pδ3jはy,zの偶関数です。さらにjのj=1,つまりxについてはこれらはゼロなのでもちろん偶関数です。

 

また,v=vy=v2,w=vz=v3はx,,zの奇関数ですから明らかに∂v2/∂xj,∂v3/∂xjはnjの偶関数です。

 

さらに,∂vj/∂x2,∂vj/∂x3についてはj=1だけの問題になりますが,uはxの偶関数ですから∂u/∂y,∂u/∂zもnjの偶関数です。以上からf2(-)=-f2(),f3(-)=-f3()が確証されました。

したがって,はx成分のみ,つまり抵抗成分のみを有すると考えられます。

 

これをDで表わすと,D=∫P1jjdS=∫PrxdS;Prx≡P1jj=P11(x/r)+P12(y/r)+P13(z/r)=-pcosθ+η[(xj/r)(∂vx/∂xj)+(xj/r)(∂vj/∂x1)] (26)となります。

そして,∂u/∂r={3aU/(4r2)}[(1+cos2θ)+a2(1-3cos2θ)/r2],∂u/∂θ={3aUcosθsinθ/(2r)}(1-a2/r2)です。

 

一方,∂/∂x=cosθ(∂/∂r)-(sinθ/r)(∂/∂θ),∂/∂y=sinθcosφ(∂/∂r)+(cosθcosφ/r)(∂/∂θ)-{sinφ/(rsinθ)}(∂/∂φ),∂/∂z=sinθsinφ(∂/∂r)+(cosθsinφ/r)(∂/∂θ)+{cosφ/(rsinθ)}(∂/∂φ)です。

 

そこで,∂u/∂x={3aUcosθ/(4r2)}[(3cos2θ-1)+a2(3-5cos2θ)/r2],∂u/∂y={3aUsinθcosφ/(4r2)}[(3cos2θ+1)+a2(1-5cos2θ)/r2],∂u/∂z={3aUsinθsinφ/(4r2)}[(3cos2θ+1)+a2(1-5cos2θ)/r2]となります。

それ故,(xj/r)(∂vx/∂xj)={3aU/(4r2)}[(1+cos2θ)+a2(1-3 cos2θ)/r2]=∂u/∂r (27)が得られます。

同様に,∂v/∂x={3aUsinθcosφ/(4r2)}[(3cos2θ-1)+a2(1-5cos2θ)/r2],∂w/∂x={3aUsinθsinφ/(4r2)}[(3cos2θ-1)+a2(1-5cos2θ)/r2] (28)より,(xj/r)(∂vj/∂x1)={3aU/(4r2)}[(3cos2θ-1)+a2(1-5cos2θ)/r2]です。

(20),(25),(27),(28)から,f1=p1jj=-pcosθ+{(3/2)ηaU/r2)}[3cos2θ+a2(1-3cos2θ)/r2] (29)です。特に,球面上(r=a)では,f1|r=a=-pcosθ+(3/2)ηU/a=(1/a){-px+(3/2)ηU}r=a (30)です。

結局,D=∫r=af1dS=2πa2-11(cosθ)f16πηaU (31)を得ました。すなわち,球は流れの向きに大きさ6πηaUの力を受けるという有名なストークスの法則が得られました。

実は,これほど苦労して偏微分方程式の境界値問題としてマニュアル通りに解かなくても,Lambや今井功氏の流体力学の教科書にもあるように,さまざまな物理的考察や手法を交えた工夫を用いることよって,より簡単に解を求めることができます。

 

しかし,当時は物理的な意味などをあまり考えることなく,なるべく式に忠実に数理的に解いても同じ答えがちゃんと得られるか?に興味があって,あえて泥臭い計算をしたのでした。

こうした偏微分方程式の境界値問題をマニュアル通りに解くという泥臭い計算は,その後確か,気象学関係で雲などによる光の散乱の散乱断面積を計算する必要があって,実行しこともありました。

 

すなわちり,ボルン,ウォルフ著の「光学の原理」にもある"光=電磁波"の古典的なミイ(Mie)散乱の散乱断面積の式を確認する際にもやったという記憶があります。

ところで,最後にストークスの法則のガリレイの相対性原理との関係を記述しておきます。

すなわち,先の"流体+小球"の系を一様な相対速度で動く並進座標系から見た場合を考えます。

 

元の座標系での位置ベクトルを,新しい系での同じ点の位置ベクトルを'で表わします。

 

便宜上t=0 で両系が一致すると仮定すると'の関係はガリレイ変換:r't (32)で与えられます。

位置に存在している流体素片の速度()はラグランジュ微分D/dtに等しいので,'(')≡D'/dt=D/dt-D(t)/dt=()-,すなわち,'(')=()- (33)です。

 

一方,小球は元の系では静止していたので,新しい系では-なる速度を有します。

また,x'j=xj-ujtより∂/∂xj(∂x'k/∂xj)∂/∂x'k (∂t/∂xj)∂/∂t=∂/∂xjです。それ故,∂2/∂xk∂xj=∂2/∂x'k∂x'jです。特に△=△'かつgrad=grad'です。

 

さらに流速のオイラー微分は,∂(,t)/∂t=∂('+t,t)/∂t-(grad')('+t,t)となります。

そこで,ナビエ・ストークスの方程式:ρ(∂/∂t)+ρ(・grad)=-gradp+η△を,新しい座標系で見ると,ρ{∂()/∂t}+ρ{()・grad'}()=-grad'p+η△'()です。

 

ところで,先に見たように'=です。また,圧力p(,t)は運動学的な量ではないのでガリレイ変換のスカラーです。すなわちp'(',t)= p(,t)です。

 

結局,ρ(∂'/∂t)+ρ('・grad')'=-grad'p'+η△''となります。

よって,新しい系でもナビエ・ストークスの方程式の形は保存されます。また,連続の方程式が保存されること:div''=0 も自明です。

したがって,ガリレイ変換(32),(33)に対してナビエ・ストークスの方程式が保存されるので,ニュートン流体に対してはガリレイの相対性原理が成立します。

また,境界条件が与えられた場合にもガリレイの相対性原理が無矛盾であることは自明です。

 

それ故,ストークス近似の流体方程式の組:η△v=gradp (1),div=0 (2) において,無限遠で静止した流体中を一様速度-=-Uで運動する小球がある系に対する流体速度の解は,既に求めた(21)~(23)式のに-を加えたもので与えられることがわかります。

以上の考察から,静止流体中を速度Uで運動する小球があるような系に対して,境界条件:=U at r=a (34),=0 ,p=p at r=∞ (35)を満足する(1),(2)の解は次式で与えられます。

すなわち,p=p-(3/2)ηaUx/r3 (36),また,u≡vx{aU/(4r)}(3+a2/r2)+{3aUx2/(4r3)}(1-a2/r2) (37),v≡vy{3aUxy/(4r3)}(1-a2/r2) (38),w≡vz{3aUxz/(4r3)}(1-a2/r2) (39)です。また,D=-6πηaU (40)も自明です。

つまり小球は進行方向に正反対の向きに大きさ6πηaUの抵抗を受けることになります。

 

ここで,レイノルズ数:R=ρUL/η=UL/νにおいて通常は代表長さLとして小球の直径d=2aを選びます。

 

このとき無次元化された抵抗係数CD≡D/(ρU2/2)=D/(ρU2πa2/2)をストークスの抵抗法則で表わすとCD24/Rとなります。

しかし,こうしたストークス近似が可能で抵抗が運動速度に比例するというストークスの抵抗法則が有効なレイノルズ数Rの範囲は非常に狭く,この法則はせいぜいR~100前後でしか成立しません。

 

通常の雨滴落下や冷却塔の温排気の凝結した水滴落下の終末速度の計算などにおいても,経験上はほとんど適用不可能であることがわかっています。

 

次は,もう少しはましな近似であるオセーン近似(Oseen近似)の計算をする予定です。

 

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

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

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

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

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

 

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

2007年7月27日 (金)

遅い粘性流(1)(Stokes近似)

学生時代は素粒子論が専門でしたが,就職難の1977年に流体の数値シミュレーションが主体の会社に入ってからは,必要に迫られて流体力学と数値計算の勉強をするようになり,これが第2の専門となりました。

 

応用物理学はやったことがなかったし,仕事もひととおりの公式類を知っていれば特に支障はなかったのですが,持ち前の好奇心もあって入社して数年間は流体力学関係の勉強をしました。

もっとも,当時はプログラマー残酷時代といわれる時代で,残業も半端ではなく,暇があると酒を飲んで発散していました。

 

それでも,寝る前などに細々と勉強して書き綴ったのが,この記事のネタ本である1980年当時の覚書きノートです。

もっとも,私は,昔から勉強が好きだったというわけではありませんでした。

 

最近,TOKIOの山口達也君が父親役をやっているドラマで,"受験の神様"と呼ばれる女子中学生に「あなたは勉強したいの,それとも受験がしたいの?」と問われてキョトンとしているシーンがありましたが,「勉強=受験」であった頃は成績は悪くなかったけど,強いられる感じの勉強は好きではありませんでした。 

試験の点数がいいとか速く解けるか?ということで他人との競争に負けるのは大して気にならなかったし,答案が不完全とかケアレスミスなども二の次でしたが,いざ自分が取り掛かった問題が解けないとか理解できないということは我慢できない性格でした。

 

試験ではそうはいかないでしょうが,わからないときは持ち帰って,ちゃんと解けて,理解できるまでは粘ったという記憶があります。 

 

未だに「高々知っているか知らないかの違いじゃないか」というような記憶するだけの知識レベルについては,あまり興味はないのですが,自然現象に限らず社会現象についても,なぜそうなるのか?というような合理性を追求することについては,理解できて解決できないと我慢できないという傾向は少なからず残っています。

 

その他,言わせてもらえば,語学のような学問は,会話で英語などの外国語とかが読めたり話せたりするだけであれば,その国のネイティブな幼児のほうがより優れていると思えますから,バイリンガルであることにこそ,その意味があると思っています。

社会人になってみると,無理に勉強しなくてもよいし,むしろ勉強や研究をする機会が少なくなりましたが,逆に好奇心だけでいろいろとやってみようという気になりました。 

余談はさておき,流体運動の現象というのはこの地球上の表面では薄い層を除いて大部分の領域では粘性が無視できるレイノルズ数(Reynolds number)の大きい流れが主体です。

 

こうした領域では完全流体の近似が成り立ちます。そして多くは乱流に移行する現象になっています。

 

しかし,ここでは,1980年当時の覚書きノートから,非圧縮かつ密度が一様な単成分流体のナビエ・ストークスの方程式(Navier-Stokes equation)において,レイノルズ数が比較的小さくストークス(Stokes)近似やオセーン(Oseen)近似が有効なケースについて流体中の球状物体に及ぼされる力=抵抗力を計算してみます。

今日はまず古典的なストークス近似の扱いからはじめます。

非圧縮,かつ密度が一様な単成分ニュートン流体のナビエ・ストークスの方程式はρ(∂/∂t)+ρ(・grad)=-gradp+η△で与えられます。

 

また,この場合の"連続の方程式=流体質量保存)の方程式"はdiv=0 です。ここでρは流体の密度,(,t)は流体速度,pは圧力,ηは第一粘性率です。

そして,レイノルズ数をR≡ρUL/η=UL/νで定義します。ν=η/ρは動粘性係数と呼ばれ,UとLは,それぞれ流れを特徴付ける代表流速と代表長さです。

 

レイノルズ数:Rが小さい場合は遅い粘性流として,ナビエ・ストークスの方程式の流速の2乗程度の"左辺第2項=慣性項"は非常に小さいので,"右辺の項=粘性項"と比較して省略する近似が可能で,これをストークス近似といいます。

特に,流速が時刻tに依らない定常流に対する定常ストークス近似の問題は,η△=gradp (1)とdiv=0 (2)を連立させて解く問題になります。

まず,(1)式の発散を取り(2)を適用すれば,△p=0 (3)を得ますからpは調和関数であることがわかります。

無限遠で,が一様速度に,圧力pがp(一定)となる非圧縮性流体中に小半径aの球を静止固定させた系において,U≡||が小さい低レイノルズ数の線型粘性流に対し定常ストークス近似を適用します。

の向きをx軸の正の向きに取り,また小球の中心を原点にとります。さらに極座標(r,θ,φ)に対してはx軸を極軸に取ります。

境界条件は次のように与えられます。0 at r=||=a (4),=Ui,p=p at r=∞ (5)です。

△p=0 (3)の解p()を(1)に代入すると,(4),(5)の境界条件に対して(1),(2)の解を求める問題は,右辺が既知となったポアソン(Poisson)方程式に対する条件付きのディリクレ問題となります。

ここで,η△G(,')=δ(') (6),G(,')=0 at r=a, ∞ (7) によって演算子η△に対する第一種のグリーン関数G(,')を定義すると,(1)の解は次式で与えられます。

 

()=∫r'=ar'=∞(,')grad'p(')d3'+ηU∫r'=∞{∂G(,')/∂r'}dS' (8)です。

なぜなら,グリーンの公式:{()△G(,')-G(,')△()}d3r=∫{()∂G(,')/∂n-G(,')∂()/∂n}dSにおいて,両辺にηを乗じη△()=gradp(),η△G(,')=δ(')を代入すれば,S上ではG(,')=0 なので(')=∫(,')gradp()d3r+η∫{∂G(,')/∂n}dSとなるからです。

ここで,積分変数を'として'を交換すれば()=∫'(',)gradp(')d3'+η∫'{∂G(',)/∂n'}dS'です。

 

さらにグリーン関数の対称性からG(',)=G(,')であり,よって明らかに∂G(',)/∂n'=∂G(,')/∂n'ですから,結局()=∫(,')grad'p(')d3r'+η∫'{∂G(,')/∂n'}dS'が得られます。

 

ただし,V'はr'=aとr'=∞で囲まれる同心球の領域はr'=aとr'=∞ の同心球面を示しており,∂/∂n'はr'=aでは内向き,r'= ∞では外向きの法線方向の微分係数を示しています。

したがって(6),(7)を満足するG(,')を求めれば,解()が得られることになります。

そこで,まず,△p=0.を満足する解p()を求めます。

 

系の極軸に対する対称性からpはrとcosθのみの球関数によって表わされます。

 

r→ ∞でp=pとなる解の一般形はp=p+Σl=0(Al/rl+1)Pl(cosθ) (9)で与えられます。ここでPl(z)はルジャンドル多項式(Legendre多項式)です。

次に,gradp=(∂p/∂x)+(∂p/∂y)+(∂p/∂z)を極座標に変換します。

 

x軸を極軸としているので,∂/∂x=cosθ(∂/∂r)+(sin2θ/r){∂/∂(cosθ)},∂/∂y=sinθcosφ(∂/∂r)-(cosθsinθcosφ/r){∂/∂(cosθ)}-{sinφ/(rsinθ)}(∂/∂φ),∂/∂z=sinθsinφ(∂/∂r)-(cosθsinθsinφ/r){∂/∂(cosθ)}-{cosφ/(rsinθ)}(∂/∂φ)です。

故に,∂p/∂x=-Σl=0{(l+1)Al/rl+2}Pl+1(cosθ) (10),∂p/∂y=-Σl=0(Al/rl+2)P1l+1(cosθ)cosφ (11) です。

 

ただし,P1l+1(cosθ)≡sinθ{dPl+1(cosθ)/d(cosθ)}={(l+1)/sinθ}{Pl(cosθ)-cosθPl+1(cosθ)}を用いています。

 

また,∂p/∂z=-Σl=0{(l+1)Al/rl+2}P1l+1(cosθ)sinφ (12) です。

 

(実際のノートには泥臭い計算の詳細があるのですが,ここでは割愛しています。)

次に(6),(7)を満足する第一種グリーン関数G(,')を求めます。

良く知られているように,△(1/|'|)=-4πδ(')です。

 

そこで,r≧aの全領域で'を除いて特異点を持たず有限確定値を取るポアソン方程式の任意のグリーン関数は,-1/(4π|'|)にラプラス方程式の一般解を加えることによって得られます。

 

ただし,この場合の極軸(対称軸)は'を結ぶ有向直線です。

 

これはグリーン関数の対称性によるものです。そこで,G(,')={-1/(4πη)}{1/|'|-Σl=0(cl/rl+1)Pl(cosω)}(r≧a)と書けます。

 

ωは'のなす角です。=(r,θ,φ),'=(r',θ',φ')なら,cosω=cosθcosθ'+sinθsinθ'cos(φ-φ'),cosω=(,')/(|||'|)が成立します。

 

また,式Pl(cosω)=Pl(cosθ)Pl(cosθ')+2Σm=1l{(l-m)!/(l+m)!}Pml(cosθ)Pml(cosθ')cos{m(φ-φ')}も成立します。

1/|'|=1/(r2+r'22rr'cosω)1/2をPl(cosω)で展開すると,r=aでG(,')=0 となるような展開係数clがPl(z)の直交性から一意的に定まります。

 

これは静電気学における鏡像の原理と同じ結果を与えます。つまりアポロニウスの定理により,G(,')の表現式において右辺第二項の級数和-Σl=0(cl/rl+1)Pl(cosω)は'と同一の向きで半径がa2/r'の位置に大きさ-a/r'の源を持つニュートン・ポテンシャルに一致します。

すなわち,G(,')={-1/(4πη)}[1/(r2+r'22rr'cosω)1/2(a/r'){1/(r2+a4/r'22ra2/r'cosω)1/2}]です。

 

そこで,G(,')={-1/(4πη)}Σl=0{rl/r'l+1-a2l+1/(rr')l+1}Pl(cosω)(a≦r<r'),G(,')={-1/(4πη)}Σl=0{r'l/rl+1-a2l+1/(rr')l+1}Pl(cosω) (a≦r'<r) (13) となります。

(13)によって,(6),(7)を満足するグリーン関数が完全に決定されたので()=∫r'=ar'=∞(,')gradp(')d3'+ηU∫r'=∞{∂G(,')/∂r'}dS' (8)に従って速度()を求めます。

 

まず,右辺第二項,すなわちvxの剰余的な項を計算します。 

以下の計算では,Pl(z)の直交性∫-11dzPl(z)Pl'(z)=2δll'/(2l+1)や,∫-11dz{Pml(z)}2{2/(2l+1)}{(l+m)!/(l-m)!},∫0dφ'cos{m(φ-φ')}=0 ,∫0dφ'cos{m(φ-φ')}cos{m'(φ-φ')}=∫0dφ'sin{m(φ-φ')}sin{m'(φ-φ')}=πδmm'等の公式を用います。

 

ηU∫r'=∞{∂G(,')/∂r'}dS'において,∂G(,')/∂r'を考えると(13)の右辺の級数展開で,この積分への寄与がゼロでないものはl=0 に対応する項のみなので次式が得られます。

 

すなわち,ηU∫r'=∞{∂G(,')/∂r'}dS'=-U・limR→∞[R2{-1/R2+a/(rR2)}]=U(1-a/r)です。

途中ですが疲れたので今日はここまでにします。

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

2007年7月26日 (木)

Priestleyの実験とCoulombの法則

 コーヒー・ブレイクとして,今日は1982年2月の私の覚書きノートから,ごく軽い話題を1つ取り上げます。

 すなわち,プリーストリーの実験(Priestleyの実験)の結果からクーロンの法則(Coulombの法則)を説明することを考えます。

問題提起としては次のようになります。

 

"一様に帯電した金属球殻の内部には電気力が存在しないとする。

 

または,球内の任意の位置に(球殻と絶縁された中空に)帯電体を入れたとき,球殻の電荷分布はそれによって影響を受けないとする。

 

こうした事実=内部の帯電体が如何なる電気力も受けないという事実,あるいは,Priestleyの実験結果だけに基づいて,

 

2つの電荷の間に働く電気力があるとすれば,その力の方向は2つの電荷を結ぶ向きであり,力の大きさはそれらの間の距離rの2乗に反比例する,という法則を証明できるか?" です。

 

まず,金属球殻についてのこの問題とは全く関係なく周りに全く何もない空間内に,2つの単位電荷1,2がそれぞれ位置ベクトル1,2で表わされる位置にあるとします。

 

このとき,2が1によって受ける力21は作用反作用の法則から,21に平行,または反平行な向きを持ち,力の大きさはそれらの間の距離|21|のみに依存すると考えられます。

 

つまり,fをある関数として21=f(|21|)(21)となるはずです。

そこで,単位電荷1を固定して,その位置を原点にとると,21にある単位電荷2が1によって受ける力=f(r)と表わすことができます。

 

このとき,力の回転を取ると,(∇×)i=εijkj(f(r)xk)=εijkjkf+(xjk)df/dr}=0 となりますから,rot=∇×=0 です。つまり,は保存力です。

したがって,ポテンシャルU(r)が存在して=-gradU(r)=-∇U(r)と書くことができます。

 

(今だったらポアンカレの補題によって"閉形式は完全形式である。"などと,少しは気のきいたことを書いたかもしれないですね。(2006年10/21の記事「ポアンカレの補題」参照))

次に,面電荷密度σを有する半径aの一様球殻がある場合の電気力を与えるポテンシャルをφ(r)とすれば,φ(r)はU(r)により,φ(r)=σ∫r'=adS'U(|'|)と表わされるはずです。dS'は球殻の面積要素です。

 

これは,謂ゆるポテンシャルの重ね合わせの原理ですね。

ここで,この球殻の球の中心を原点とし極角θを'のに対する偏角とすれば,φ(r)=2πa2σ∫-11(cosθ)U({r2+a22racosθ}1/2)=(2πaσ/r)∫-r+ar+aRU(R)dRとなります。

 

最初から,rだけの関数であろうと予想してポテンシャルをφ(r)と書いた通り,右辺はrのみの関数になっています。

そして,"仮定=Priestleyの実験結果"によれば, 0=-gradφ=-∇φ=-(dφ/dr)(/r)(r<a),よってdφ/dr=0 です。

 

したがって,d2(rφ)/dr2=d(rdφ/dr+φ)/dr=dφ/dr=0 となることが必要であるとわかりました。

そこで,V(R)≡RU(R)と置くとrφ(r)=2πaσ∫-r+ar+a(R)dRより,V'(r+a)-V'(-r+a)=0 が成立します。ただしV'(r)≡dV/drです。

 

よって,dV/dr=(定数)です。そこでC0,C1を積分定数としてV(r)=rU(r)=C0r+C1,つまりU(r)=C0+C1/rと書くことができます。

 

V(r)=rU(r)=C0r+C1を,rφ(r)=2πaσ∫-r+ar+a(R)dRに代入すると,φ(r)=(定数)になりますから確かに恒等的にdφ/dr=0 です。

 

したがって,()=-gradU(r)=-∇U(r)=(C1/r2)・(/r)∝ /r2 (ただし/r)が得られました。

 

すなわち,電気力()の方向はと同じで大きさは距離rの2乗に反比例します。

 

ただし,()が引力か斥力か?は,これだけからはまだ不明です。

 

PS:クーロンの法則を実験的に発見したのは名称の通りクーロン(Coulomb)であるとされています。

 

しかし,プリーストリ-(Priestley)が"静電気力は距離の二乗に反比例する"という性質を予測して本記事のような実験をしたのは彼よりも前のことです。

 

また,プリーストリーより後ですがクーロンより前に独立に精密な定量的実験をしてこの法則を検証したのはキャベンディッシュ(Chavendish)でした。

 

時代が時代ならプリーストリー・キャベンディッシュ(Priestry-Chavendish)の法則となっていたかも知れませんね。

 

キャベンディッシュ(Chavendish)は万有引力(重力)の法則においても逆二乗則が成り立つことを実験的に検証しました。

 

実験の詳細については,例えば「FNの高校物理(キャベンディッシュの重力実験)」を参照してください。

 

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

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

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

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

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

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

2007年7月25日 (水)

中性子の磁気モーメント

今日は昔の学生時代を思い出して少し素粒子論を復習してみよう

と思います。

 
現在ではハドロン(Hadron)を構成するクォーク(Quark)の

フレーバー(Flavor)の自由度は,2種類ずつ3世代あって,

アップ(up),ダウン(down);チャーム(charm),ストレンジ

(strange);ボトム(bottom),トップ(top)の6種類があると

されていますが,

 "核子=陽子p or 中性子n"を構成するだけなら,

"アップ=u"と"ダウン=d"の2種類だけで十分です。

 

 核子と,その核力を媒介する主要な粒子であるπ中間子を

uとd,あるいはこれらの反粒子の複合粒子として記述する

には,2次元特殊ユニタリ群:SU(2),すなわち,アイソスピン

(Isotopic spin;荷電スピン)の群の表現を考えるだけで

十分です。


 
"核子=陽子p or 中性子n"は,uとdの3体で構成される

バリオン(Baryon;(重粒子)であり,スピンが1/2のFermi粒子

(Fermion)です。

 
そして素電荷eを電荷の単位とすると,陽子pは電荷が1

であり,中性子nは,電荷がゼロである,ということによって

特徴付けられています。

 そして電荷をQ,アイソスピンベクトルを,バリオン数をBと

すると,Q=I3+B/2という等式が成立します。

 核子はそのアイソスピンがI=1/2で,陽子はI31/2,中性子

はI3=-1/2の固有状態です。

 
ただし,ストレンジネス(Strangeness):Sまで考慮に入れると,

Q=I3+Y/2;Y≡S+Bとなります。

 Yはハイパーチャージ(Hyperchage)といわれる量です。

 私くらいの世代の学生時代なら,丁度卒業の頃くらいにチャーム

(Charm)を与える(J/ψ)粒子が発見されたわけですから,わかって

いたのはこの程度までです。

 
ちなみにスピンがゼロの擬スカラー粒子であるπ中間子に

 ついては,B=0,I=1であり,π±0は,それぞれ,

 I3=±1,0 の固有状態です。

 一般に各素粒子はユニタリ群の既約表現の1つ1つに対応

 していて,超選択則(Superselection rule)により,崩壊現象

 を除けば異なる既約表現間の遷移は禁止されていると考え

 ます。

 一方,クォークu,dのスピンは1/2で,それらの電荷Qは,

 eを単位として,それぞれ,2/3,-1/3です。

 ということは,B=1/3なので,Q=I3+B/2からu,d

 のアイソスピンはI=1/2で,それぞれはI3がI31/2,

 -1/2の固有状態であるとしてよいことがわかります。

 
SU(2)群の2次元基本表現を表わすクォーク3体で構成される

 複合粒子=核子について,3体が合成された直積表現の既約表現へ

 の分解は, 2×2×2=2+2+4 となります。

(因みにストレンジネス:Sまで含めたSU(3)だと,

3×3×3=1+8+8+10 になります。)

 

さて,まず,2体での既約分解が,2×2=1+3となるのは自明です。

 

つまり,Tij=T[i,j]+T{i,j};T[i,j]1/2(ij-Tji)(反対称=1)

{i,j}1/2(ij+Tji)(対称=3) ですね。

 
そして,2×2=1+3 から,さらに3体では,

2×2×2=(1+3)×2=2+(2+4)という既約分解

を得ます。

 
まず,右辺の最初の2はアイソスピン 0 と 1/2を合成して

アイソスピン 1/2を作ることに相当します。


 次の2はアイソスピン1と1/2を合成してアイソスピン1/2を

つくることです。

 

右辺の最後の4はアイソスピン 1と1/2を合成してアイソスピン

3/2をつくることに相当します。

 
アイソスピン3/2はTijkを完全対称にすることで得られますから

[i,j,k]です。

 このテンソルの成分の数はi,j,k=1,2の全ての組み合わせ

の数,つまり,2個から重複を許して3個を選ぶ組み合わせ

なので確かに4個あります。

 
これは核子ではなくて,π-pあるいはπ-n共鳴であり,

質量が,およそ,1230MeVのΔ粒子(デルタ)を表わしています。

一方,2×2×2=2+(2+4)の右辺の最初の2はT[I,j]k,

すなわち,(1/2)(121-T211)と(1/2)(122-T212)を表わして

います。

よって右辺の残りの2はT{I,j}K[I,j,k]で与えられます。

これらのゼロでない独立な成分は,(1/3)(2112-T121-T211)

(1/3)(2221-T212-T122)の2つです。

そして,陽子pと中性子nは,この最後に示した方のI31/2,

-1/2の既約表現に対応するとされています。

 
これらはu,dという記号をそのままu,dを示す状態の波動関数

として表現すれば,規格化も含めて,

p=(1/6)1/2(2uud-udu-duu),

n=(1/6)1/2(2ddu-dud-udd)

となります。

 
ところで,こうした理論によるとアイソスピンとスピンが共に3/2

Δ++ 粒子において,スピン成分がsz=+3/2の状態は,

Δ++=u↑ となります。

 

これはフレーバー自由度,スピン自由度について共に

完全対称です。

 
ハドロンのクォークによる複合粒子としての表現が"軌道角運動量

がゼロの基底状態=S状態"で与えられるという仮定によれば,

Δ++はクォークの交換に対して位置座標の交換を含めて完全対称

な状態関数で表現されることになります。

 
しかし,これはFermi統計,つまり多粒子系の状態はFermi粒子の

交換に対して反対称であるべきである,という要請に矛盾します。

そこで,実際の理論では,もう1つ別の自由度であるカラー(Color)

というものを導入し,カラー自由度については1重項(無色:Singlet)

であること,つまりクォークの交換についてカラー自由度について

完全反対称の状態にあるとして,この矛盾を解消しています。

そこで,今問題としている陽子:

p=(1/6)1/2(2uud-udu-duu)と,

中性子:n=(1/6)1/2(2ddu-dud-udd)

について考えると,これらは1番目と2番目のクォークの交換

について対称です。

 
しかし,カラー自由度については完全反対称ですから,スピン

の自由度についても1番目と2番目のクォークの交換について

対称であることが要求されます。

そこでスピン1/2に対する回転群SU(2)の既約表現についても

同じ変換性を持つ表現で,

|↑>=(1/6)1/2(2↑↑↓-↑↓↑-↓↑↑)と,

|↓>=(1/6)1/2(2↓↓↑-↓↑↓-↑↓↓)

を採用します。

陽子のスピンアップ状態としては,

(1/6)(2uud-udu-duu)(2↑↑↓-↑↓↑-↓↑↑)

と表現すればいいのでは?と推測されます。

そこで,結局アイソスピンとスピンの両方を考慮したときに,

クォークの交換に対して完全対称でなければならないことから

|p>=(1/18)1/2[uud(2↑↑↓-↑↓↑-↓↑↑)+cyclic]

=(1/18)1/2[2u-u-u)+cyclic]

と表現さるべきであることが結論されます。

 
同様に,

|n>=(1/18)1/2[ddu(2↑↑↓-↑↓↑-↓↑↑)+cyclic]

=(1/18)1/2[2d-d-d)+cyclic]

です。

 
これらはかつて流行したことのあるフレーバー・スピン対称性;

SU(6)の対称な,56重項既約表現に対応するものです。

余談ですが,このSU(6)というのは,対称性の帰結として

角運動量保存則が従う現実世界が等方的であるという性質,

つまり現実の空間での回転群:SU(2)~SO(3)というスピン

角運動量に関わる対称性と,フレーバーという内部空間,

アイソスピンならアイソ空間(荷電空間)での回転対称性

を含むSU(3)のフレーバー対称性を合成したものです。

こうした現実空間と仮想内部空間を混合した対称性というのは,

非相対論で成り立つ近似的なものです。

 
こうした混合対称性は相対論まで含めた4次元時空という実空間

に対しては厳密には成立し得ないことが定理として証明されて

います。

,

ただし,例外があって超対称性はこの限りではありません。
(※ワインバーグの「場の量子論」邦訳第5巻参照※)

クォークは電子などのレプトン(Lepton;軽粒子)と同じく,

構造を持たないスピンが1/2の素粒子なので,その磁気回転比

gは,g=2で近似することができます。

そこで,クォークで構成された複合粒子の磁気モーメント

(磁気能率)をμとし,電荷を持つ構成粒子によって,これを評価

すると,z成分は,

μzΣi{eic/(2Mi)}(liz+giz)

で与えられると考えられます。

 

ここで,ei,Mi,i,iは,それぞれi番目の構成粒子の電荷,質量

軌道角運動量,スピンです。

 また,h
c≡h/(2π)はPlanck定数,cは光速です。

iz=σi3/2と書き,M~ Muとすると,

μz|p>=ehc/(2Mu)(1/18)1/2

[{(10/3)u(1/3)u(1/3)u}

+{(10/3)u(1/3)u(1/3)u}

+..]となります。

したがって,普通に期待値として陽子の磁気モーメント

μpz=<pz|p>を計算すれば,
μ
pz{ehc/(2Mu)}×3×(1/18)×[(5/3)×4+(-1/3)+(-1/3)]

=ehc/(2Mu) が得られます。

 

同様に中性子では,

μnz=<nz|n>={ehc/(2Mu)}×3×(1/18)

×[(-4/3)×4+(2/3)+(2/3)]=(-2/3){ehc/(2Mu)}

が得られます。
 

したがって,理論的には,中性子と陽子の磁気モーメントの比

として,μn(-2/3)という結果を得ます。

実験によると,Bohr磁子μB=ehc/(2mp)を単位として,核子

の磁気モーメントは,陽子pがおよそ2.79で,中性子nが-1.91

であることがわかっています。

 
つまり,μp 2.79μBn ~ -1.91μBです。

 
実測値では,中性子と陽子の磁気モーメントの比は,

μnp ~ -1.91/2.79 ~ (-2/3)で与えられること

になります。

それ故,先の理論的考察は実験事実を正しく評価しています。


 
普通,電荷を持たない物体では角運動量があっても,電流がない

ので磁気モーメントはゼロですから,これは中性子が総体として

の電荷はゼロでも,内部に電荷密度の分布を持つ構造がある証拠

を与えていると考えられます。

さらに,μp がBohr磁子μB2.79倍であることが,陽子の質量mp

がクォークの質量Mu2.79倍程度であることを示唆している

と考えるなら,up,downクォークの質量がMu ~ Md 330MeV

であろうというクォーク質量の推定値も得られます。

 参考となる書籍が,徒歩では行けないちょっと離れたトランク

ルームにあり,おまけに風邪で外出もままならないので,ほとんど

記憶に頼って計算したので,合理的な結果を得るまでかなり計算

間違いを繰り返して苦労しました。

参考文献:S.Weinberg著(青山秀明他共訳)「場の量子論5」

(超対称性;構成と超対称標準模型)(吉岡書店)

ランキングへ ← クリックして投票してください。

(1クリック=1投票です。1人1日1投票しかできません。)

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

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

(ブログ村科学ブログランキング)

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

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

2007年7月23日 (月)

e とπの超越性

 木曜からの風邪がだんだんひどくなり,セキが止まらないので仕方なく近くの滝野川病院で診察を受け,薬をもらいました。

 近くなのに病院まで歩いていけなくてタクシーを拾いました。普通は風邪くらいなら市販薬を飲んで寝るだけなのですが,手術の予後で体力回復が万全でないので心配になったのです。

 それにしても私はこのところ病院づいていますね。

しかし,もらった薬は食後に飲むように指定されていて,自宅には米以外の食料がなかったので病院のすぐそばにあったスーパーでいろいろと食料を買いましたが,自宅まで歩けないのでまた仕方なくタクシーに乗りました。

 

でも自宅に帰ってもおかゆさえ食べる気にならず,まあいいやと食前だけど薬を飲んだところ,プラシーボかもしれないが,さすがは病院の薬であるせいなのか,すぐに症状は軽減されました。

最近,ブログネタもないので,本の受け売りですが,以前の2007年5月27日の記事「代数的数と超越数」の続きとして,eとπの超越性の証明でも書いておきます。

(eの超越性の証明)

eが代数的数であると仮定すると,a0+a1e+..+ann0 を満たす整数の組a0,a1,..,anが存在します。n≧1,a0n0 です。

ところで,f(x)をm次の多項式,tを複素数とし積分路を0とtを結ぶ直線としてtの関数I(t)をI(t)≡et0t-x(x)dxで定義し,部分積分を繰り返すと,I(t)=et(0)-f(t)+et0t-x(df/dx)dx=etΣj=0m(j)(0)-Σj=0m(j)(t)が得られます。

(x)をf(x)≡xp-1(x-1)p..(x-n)pとします。

 

pは十分大きい素数です。

 

J≡-{a0(0)+a1(1)+..+an(n)}とおくと,J=-(a0+a1e+..+annj=0m(j)(0)+a0Σj=0m(j)(0)+a1Σj=0m(j)(1)+..+anΣj=0m(j)(n) =a0Σj=0m(j)(0)+a1Σj=0m(j)(1)+..+anΣj=0m(j)(n)となります。

ただし,m=degf=np+p-1であり,0≦j≦mなるjに対してf(j)(x)はxp-1,(x-1)p..,(x-n)pの各々をそれぞれ0,1,..,n回ずつ(Σi=0li=j)微分した項の積和ですから,右辺の各f(j)(k)(k=0,1,2,..,n)は全て整数であり,それゆえJも整数です。

 

そして,j<pなら明らかにf(j)(k)=0 (k=1,2,..,n)でj<p-1ならf(j)(0)=0 です。

ところで,k≧lのときk(k-1)..(k-l+1)/l!はl!で割り切れるので整係数多項式g(x)のl回導関数の係数はすべてl!で割り切れます。したがってj≧pのときf(j)(x)の係数はj!で割り切れ,よって全てp!で割り切れます。

したがって,f(j)(k)はj=p-1かつk=0 の場合を除いて,p!で割り切れます。ところが,f(p-1)(0)=(-1)p(p-1)!(n!)pであり,これは(p-1)!では割り切れますが,pが十分大きいときpでは割り切れません。

 

それ故,Jはゼロでない整数であり,(p-1)!で割り切れます。よって,|J|≧(p-1)!です。

一方,0≦x≦k,0≦k≦nのとき,|x|,|x-1|,..,|x-n|≦2nなので |I(k)|=|∫0kk-x(x)dx|≦kek(2n)です。

 

それ故,|J|≦|a0||I(0)|+|a1||I(1)|+..+|an||I(n)|≦(2n)Σk=0n|ak|kekです。m=np+p-1ですから,pに無関係な定数c>0 が存在し,|J|<cp が成立します。

ところが,Stirlingによれば,pが大きいとき,(p-1)!~pp-pですから,素数pを大きくすると,やがて(p-1)!>cpとなります。

 

そこで,2つの不等式|J|≧(p-1)!と|J|<cp は両立しませんから,これは矛盾です。したがって,eが代数的数であるとした仮定が間違いであり,eは超越数であることがわかります。(証明終わり)

(πの超越性の証明):

πを代数的数であると仮定すると,iπも代数的数です。ここでθ=iπの共役:つまりθ=iπの最小多項式における相異なるn個の根の組をθ1=θ,θ2,..,θd;d≡degθとします。

 

このとき等式(Euler's formula)e1=0 より,(θ11)(θ21)..(θd1)=0 となります。

 

そして左辺を展開すると,Σi=1,2,..,nδi=0,1δ1θ1+δ2θ2+..+δdθd0 となります。

2d個のδ1θ1+δ2θ2..+δdθdi0または1)のうちで 0 でないものがn個あったとし,それらをα1,α2,..,αnをとします。

 

また,n<j≦2d に対してはαj0 とします。そうすると,q+α1α2+..+αn0 が得られます。ここで,q2d0です。

ここで,eの超越性の証明と同様に,I(t)=et0t-x(x)dxを考え,今度はf(x)≡lnpp-1(x-α1)p..(x-αn)pとします。

 

ただしpは十分大きい素数で,lはlθ1,lθ2,..,lθd ,したがってlα1,lα2,..,lαn 代数的整数となるような整数です。

 

ここで,βが代数的整数であるとは,代数的数βの整数係数の最小多項式がその最高次の係数=1の形つまり整数係数のモニックになることをいいます。

そして,J≡-{I(α1)+I(α2)+..+I(αn)}とおくとJ=-(α1α2+..+αnj=0m(j)(0)+Σj=0m(j)(α1)+Σj=0m(j)(α2)+..+Σj=0m(j)(αn)=qΣj=0m(j)(0)+Σj=0m(j)(α1)+Σj=0m(j)(α2)+..+Σj=0m(j)(αn), m=np+p-1です。

 

多項式f(x)の各係数はlα1,lα2,..,lαnの整数係数の対称式,したがって,lθ1,lθ2,..,lθdの整数係数の対称式です。

ところが,lα1,lα2,..,lαnを根とするモニックの最小多項式をP(x)とすると,P(x)=(x-lα1)(x-lα2)..(x-lαn)=xn-s1n-1+s2n-2-..+(-1)snであって,s1,s2,..,snlα1,lα2,..,lαnの基本対称式であり,lα1,lα2,..,lαnが代数的数なのでこれらは整数です。

 

しかもlα1,lα2,..,lαnの全ての整数係数の対称式は基本対称式の整数係数の多項式として表わせることは明らかなので,(x)の各係数は整数になります。

したがって,f(j)(0)も整数です。また,Σk=1n(j)(αk)(j≧0)もlα1,lα2,..,lαnの整数係数の対称式なので,上と同じ論旨から,これは整数になります。

 

そして,j<pなら明らかにf(j)(αk)=0 (k=1,2,..,n)で,j<p-1ならf(j)(0)=0 です。したがってj≧pのとき,f(j)(x)の係数は全てp!で割り切れます。

以上から,f(j)(αk)とj≠p-1のときのf(j)(0)は,p!で割り切れます。ところが,f(p-1)(0)=(-1)p(p-1)!(lα1lα2..,l)pであり,これは,(p-1)!では割り切れますが,pが十分大きいときpでは割り切れません。

 

それ故,Jはゼロでない整数であって,(p-1)!で割り切れます。よって,|J|≧(p-1)!です。

一方,0≦x≦αkまたは,αk≦x≦0 for 0≦k≦nなるxに対して,M=maxk|2αk|とおくと,|x|,|x-α1|,..,|x-αn|≦Mなので,|I(αk)|=|∫0αkαk-x(x)dx|≦|αk|αknpです。

 

それ故,|J|≦|I(0)|+|I(α1)|+..+|I(αn)|≦lnpΣk=0n|αk|αkと評価できます。ここで,m=np+p-1ですからpに無関係な定数c>0 が存在して|J|<cpが成立します。

pが大きいとき,(p-1)!>cpとなりますから,2つの不等式|J|≧(p-1)!と|J|<cp は両立しません。したがってπが代数的数であるとした仮定が間違いでありπは超越数であることがわかります。(証明終わり)

これで,eおよびπは超越数であることが証明されました。超越数ですから当然無理数ですが,いささか本末転倒ながらこれらが無理数であることをもっと直接的に証明してみます。

(eが無理数であることの証明):

e=1+1/1!+1/2!+1/3!+..ですから,0<n!e-n!Σk=0n(1/k!) =n!Σk=n+1(1/k!)≦{1/(n+1)}(1+1/2+1/22+..)=2/(n+2)→ 0 as n→∞ が成立します。

ところが,一般にαが有理数でα=a/b(a,b∈Z,b>0)なら,αと異なる全ての有理数p/qに対して|qα-p|=|qa-pb|/b≧1/bですから,αにのみ依存する正の定数cが存在して,|α-p/q|≧c/qとなります。それ故,0<|qnα-pn|→ 0 as n→ ∞ なる整数列{pn},{qn}が存在するならαは無理数です。

以上から,eは無理数であることが示されました。(証明終わり)

(πが無理数であることの証明):

実係数の多項式f(x)に対してF(x)≡f(x)-f(2)(x)+f(4)(x)-f(6)(x)+..とおきます。f(x)は多項式なので右辺は有限項で終わりF(x)も多項式になります。

このとき,{F'(x)sinx-F(x)cosx}'={F"(x) +F(x)}sinx=f(x)sinxなので,∫0π(x)sinxdx=F(0)+F(π)と書くことができます。

今,πが有理数であると仮定して,π=a/b(a,b∈Z,b>0)とし,上式でf(x)=(bn/n!)xn(π-x)n(1/n!)xn(a-bx)nとおきます。

 

ただし,nは十分大きい正の整数とします。明らかに,j<nのとき(j)(0)=0 で,j≧nのときは多項式:{xn(a-bx)n}(j)のすべての係数は,j!したがってn!で割り切れます。故に,j≧0なるすべてのjについてf(j)(0)は整数です。

ところで,j≧0なるすべてのjについてf(j)(x)=(-1)j(j)(π-x)よりf(j)(π)=(-1)j(j)(0)も整数です。

 

以上から,∫0π(x)sinxdx=F(0)+F(π)の右辺{F(0)+F(π)}は整数です。

 

一方,左辺は,0<∫0π(x)sinxdx≦(bnπ2n/n!)∫0πsinxdx=2bnπ2n/n!<cn/n!です。cはnに依存しない正の定数です。

 

これからn! ≦cn ですが,nが十分大きいときには成立しません。

以上からπは無理数であることが示されました。(証明終わり)

参考文献:塩川宇賢 著「無理数と超越数」(森北出版)

 

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

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

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

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

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

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

2007年7月20日 (金)

揺動散逸定理

 非平衡統計熱力学の1過程としてブラウン運動などに関わる揺動散逸定理(fluctuation dissipation theorem)について述べてみます。

 

 この定理は誰が起源なのかよく知らないのですが,日本では線型応答理論の一環として統計物理学の重鎮であった久保亮五先生などが関わっていたと記憶しています。

 

 一般に,ある物理量(示量変数の組):=(a1,a2,..,an)があって,エントロピーSがの関数であるとき,系の時間発展はに対する1階微分方程式で表現されて,それはdai/dt=(dai/dt)rev+(dai/dt)irrのように,可逆な部分(dai/dt)revと不可逆な部分(dai/dt)irrの和として与えられます。

 そして,可逆部分はさらに,(dai/dt)rev=Σj{ai,aj}(∂S/∂aj)という構造を持つとします。

 

 ここで,{X,Y}はポアソン括弧のように,{Y,X}=-{X,Y}という反対称性を持ち,{X,f}=Σj{X,aj}(∂f/∂aj)という性質で規定される量であるとします。

こう定義すると,(dS/dt)rev=Σi(∂S/∂ai)(dai/dt)rev=ΣiΣj(∂S/∂ai){ai,aj}(∂S/∂aj)={S,S}=0 となって,可逆過程ではエントロピーは生成されないことになります。

 

実際には,断熱可逆変化か,あるいは可逆であって,かつサイクルである場合以外なら,可逆過程でもエントロピーの変化はありますから,これはどう解釈すべきなんでしょうか?

 

dai/dt=(dai/dt)rev+(dai/dt)irrは"可逆な部分と不可逆な部分の和である。"というよりもむしろ,"エントロピー非生成部分とエントロピー生成部分の和である。"と事実のままを述べた方がいいのかもしれません。

一方,不可逆部分については現象論的に(dai/dt)irr=Σjij(∂S/∂aj)と表わすことにします。

 

というのは,ajが示量変数のとき∂S/∂ajは示強パラメータであり,平衡の近傍では不可逆部分は示強パラメータ,あるいはその空間勾配に比例して進行するからです。

実際,問題としている系を局所平衡状態にある部分系の集まりと考えると,それぞれの部分系の物質密度をρ,単位質量当たりのエントロピーをsとしたとき,エントロピー密度(ρs)の変化は,一般に熱力学の関係式により示強パラメータFiと示量変数aiによって,d(ρs)=Σiidaiと表わされます。

 

この表式では確かに示強パラメータは,Fi=(∂S/∂aj)/Vと表わされています。

そして今,対象としている系が隣り合う2つの部分系A,Bだけから成るとし,A,Bが等しい体積Vを持つとするとき,各部分系のエントロピーS=Vρsは,Xi=Vaiの関数であると考えられます。

 

今,A,Bのエントロピーを,それぞれSA,SBとし,XiがAからBにΔXi=VΔaiだけ移動するとします。

iが系全体では保存する量であって,初めAにはXiがXiAだけBにはXiBだけあったとすると,AからBへのΔXi=VΔaiの移動による系全体のエントロピーSの増分は,ΔS=SA(XiA-ΔXi)+SB(XiB+ΔXi)-[SA(XiA)+SB(XiB)]~Σi(-∂S/∂XiA+∂S/∂XiB)ΔXi=Σi(-FAi+FBi)ΔXiとなります。

 

ここでFAi,FBiはそれぞれ部分系A,BにおけるFiの値を表わしています。

系全体を孤立系と考えると熱力学第二法則によって,ΔS>0 でなければならないので,1つの示量変数Xi=Vaiのみに着目してAからBへと微小量ΔXi>0 の移動が起こるためには,示強パラメータFiについてはFBi>FAiであることが必要になります。

このことから,示量変数Xi=Vaiが保存量のときはXiの輸送を引き起こす駆動力となるのは示強パラメータFiの空間勾配であると考えられます。

また,示量変数Xi=Vaiが非保存量のときはΔS=ΣiiΔXiにおいてFi=(∂S/∂ai)/V>0 ならば,ΔXi=VΔai;Δai>0 なる変化が不可逆過程として進行し得ます。

 

つまり,一般にdS=Σiidai,Fi=(∂S/∂ai)と書けますが,量akが保存量のとき,その保存方程式は∂ak/∂t+∇k=0 であり,その流れkは一般にk=Σjkj∇Fjと,示強的な量Fj=(∂S/∂aj)の勾配を駆動力とする形に表わされます。

 

そこで,k'≡akとおくと,k=ak=ak(d/dt)ですから,dk'/dt=d(ak)/dtk=Σjkjj=Σjkj∇(∂S/∂aj)です。

 

そこで,もしもSがajの1次関数なら,∇(∂S/∂aj)~∂S/∂j'となり,示量変数k'の求める時間発展の形式dk'/dt=Σjkj(∂S/∂j')が得られます。しかし正直なところかなり苦しいです。

こうして,はっきりと証明されたわけではないのですが,現象論的発展方程式はdai/dt={ai,S}+Σjij(∂S/∂aj)と表わされるとします。

 

これは非平衡な初期条件から平衡状態への緩和を表わしています。

 

そこで(t)の時間発展は初期条件(0)=に対してai(t)=ai+t[{ai,S}+Σjij(∂S/∂aj)]+..で与えられます。

平衡状態においても,一般に巨視的な物理量(a1,a2,..,an)はゆらいでいます。たまたま,ある時刻に(t)が(t0)=という値をとったとして,その後の時間発展を観測します。

 

(t0+t)は初期時刻t0によってさまざまな値を取りますが,それらの平均を取ったものも現象論的発展と同じになるとします。

 

すなわち,初期条件(t0)=が与えられると短い時間では平均的に<ai(t0+t)>(t0)==ai+t[{ai,S}+Σjij(∂S/∂aj)]となると仮定します。これを線型減衰の仮定と呼びます。

平衡状態におけるゆらぎの時間相関関数<ai(t0+t)ak(t0)>を求めるには<ai(t0+t)>(t0)=に初期値ak(t0)=akを掛けてakの分布について平均すればいいので,t≧0 に対して<ai(t0+t)ak(t0)>eq=<<ai(t0+t)>(t0)=keq =<aikeq+t[<{ai,S}akeq+Σjij<(∂S/∂aj)akeq]となります。

平衡状態におけるゆらぎに対するボルツマン・アインシュタインの原理,つまり,ボルツマンの原理S=kBlnWから,逆に微視的状態数WがW()=exp[S()/kB]と書ける,ことを用います。

 

全状態数をWとしたときにW()/Wが変数の状態が実現する確率となりますから,平衡状態での関数f()の平均値は<f()>eq=(1/W)∫da1..danf()exp[S()/kB]で与えられます。

 

そこで,<(∂S/∂aj)akeq=∫da1..dan(∂S/∂aj)akexp[S()/kB]=-kBδkjとなります。

したがって,<ai(t0+t)ak(t0)>eq=<aikeq+t[<{ai,S}akeq-kBik]となりますが,右辺の{ai,S}はdai/dtの可逆部分です。

 

ところが,平衡状態では物理量の任意の関数は可逆変化では変化しないので,<{f(),S}>eq=0 です。もっともこれはその現象論の範囲で厳密に証明できるわけではないので仮説として導入するわけです。

 

そして,この式でf()=aikを代入すると,<aj{aiδjk,S}+aj{aiδijk,S}>eq=0 :すなわち<{ai,S}akeq=-<ai{ak,S}>eqが得られます。

 

そこで、物理変数ai,akの時間反転対称性に関して次の2つの場合を考えます。:

(Ⅰ)変数ai,akが共に時間反転に対して対称である場合

このときは<ai(t0+t)ak(t0)>eq=<ai(t0-t)ak(t0)>eqです。さらに平衡状態の定常性から物理量の時間相関は時間差のみに依存しt0には依存しないので<ai(t0-t)ak(t0)>eq=<ai(t0)ak(t0+t)>eqです。

 

以上から,<ai(t0+t)ak(t0)>eq=<ai(t0)ak(t0+t)>eqが得られます。

そこで両辺に線型減衰の仮定を適用すると,<aikeq+t[<{ai,S}akeq-kBik]=<aikeq+[<{ak,S}aieq-kBki]となるはずです。

 

これが実際に成り立つためには,<{ai,S}akeq=<{ai,S}akeqかつLik=Lkiが満たされなければなりません。

前者は<{ai,S}akeq=-<ai{ak,S}>eqと組み合わせると<{ai,S}akeq=<{ai,S}akeq=0 となります。要するに,ai,akが共に時間反転に対して対称ならば,その時間微分{ak,S},{ai,S}は明らかに時間反転に対して反対称となりますから,時間反転対称なaiやakとの積は平衡状態ではゼロとなることを表現しています。

(Ⅱ) 時間反転に対して変数aiが反対称でakが対称の場合

このときは<ai(t0+t)ak(t0)>eq=-<ai(t0-t)ak(t0)>eqです。そこで(Ⅰ)と同様に考えて<ai(t0+t)ak(t0)>eq=-<ai(t0)ak(t0+t)>eqとなります。

 

したがって線型減衰の式から<aikeq+t[<{ai,S}akeq-kBik]=-<aikeq-t[<{ak,S}aieq-kBki](t≧0)ですが,これが成り立つためには<aik>=0 でかつLik=-Lkiが満たされなければなりません。

かくして,変数ai,akが同じ時間対称性を持つときにはLik=Lki, 反対の時間対称性を持つときにはLik=-Lkiとなります。

 

さらに係数Likが時間反転に対して反対称な外部パラメータ(は例えば磁場や速度)に依存するときには,それぞれLik()=Lki(-),Lik()=-Lki(-)となります。

発展方程式の不可逆部分を熱力学的力∂S/∂で表現する輸送係数Lijについての上述の対称性をオンサーガーの相反定理と呼び,この係数Likをオンサーガー係数と呼びます。

さらに発展方程式の不可逆部分はエントロピー生成をするという熱力学第二法則の要請から係数行列(Lij)は正値行列です。

  

なぜなら,dS/dt=(dS/dt)irr=Σi(∂S/∂ai)(dai/dt)irr=ΣiΣj(∂S/∂ai)Lij(∂S/∂aj)>0 となるべきことが要求されますが,(∂S/∂ai)はベクトルとして任意の値を取ると考えてよいからです。

平衡状態での巨視的変数の現象論的な発展方程式がdai/dt={ai,S}+Σjij(∂S/∂aj)で与えられるので,必然的に存在する"ゆらぎ=揺動あるいは雑音"Ri(t)の存在を考慮すると,の時間発展は一般的な確率微分方程式dai/dt={ai,S}+Σjij(∂S/∂aj)+Ri(t)で表わされると考えられます。

 

ここで通常はゆらぎRi(t)は完全にランダムであり<Ri(t)>=0 ,<Ri(t)Rj(t')>eq=2Dijδ(t-t')なる白色雑音(white noise)で与えられます。

こうすると,時刻tにおいて状態が実現する確率P(,t)に対して,次のフォッカー・プランク方程式(Fokker-Planck)が得られます。

  

∂P(,t)/∂t=-Σi(∂/∂ai){ai,S}P(,t)+ΣiΣj(∂/∂ai)[-Lij(∂S/∂aj)+Dij(∂/∂aj)]P(,t)です。

一般に,1次元で考えたとき外力Fがないときのゆらぎも含めた粒子の運動は,その速度をuとするとき,次の"運動方程式=ランジュバン方程式(Langevin)"du/dt=-γu+R(t)/m に従います。

 

そして,ここでもゆらぎR(t)は白色雑音,つまり<R(t)>=0 ,<R(t)R(t')>=2Duδ(t-t')を満たしているとします。

このとき,時刻t1に速度u1を持っていた粒子が時刻tに速度uを持つ条件付の確率分布T(u,t|u1,t1)は,(∂/∂t)T(u,t|u1,t1)=γ(∂/∂u)[u+(kBT/m)∂/∂u+(Du/m2)∂2/∂u2]T(u,t|u1,t1)という方程式に従うことがわかります。

 

確率分布が従うこの方程式を,フォッカー・プランク方程式と呼ぶわけですね。

そして,先のについての運動方程式dai/dt={ai,S}+Σjij(∂S/∂aj)+Ri(t)を上の1次元速度uに対するランジュバン方程式におきかえ,時刻tに状態が実現する確率分布P(,t)を上の確率分布T(u,t|u1,t1)におきかえれば,先述のP(,t)に対するフォッカー・プランク方程式が得られます。

これが定常解として平衡分布Peq(a)=exp[S()/kB](あるいはこの定数倍)を持つためにはkBij=Dijとなることが必要条件になります。

つまり,∂Peq()/∂aj=(1/kB)(∂S/∂aj)Peq()なので∂P(,t)/∂t=-Σi(∂/∂ai){ai,S}P(,t)+ΣiΣj(∂/∂ai)[-Lij(∂S/∂aj)+Dij(∂/∂aj)]P(,t)においてP(,t)=Peq()とおくと,kBij=Dijが満たされる場合には右辺の第2項はゼロになります。

一方,第1項はΣi(∂/∂ai){ai,S}Peq()=Peq()[Σi(∂/∂ai){ai,S}+(1/kBi(∂S/∂ai){ai,S}]となりますが,定義によってΣi(∂S/∂ai){ai,S}={S,S}=0 です。

 

これほど自明ではありませんが,Σi(∂/∂ai){ai,S}=0 も成立します。

実際,Σi(∂/∂ai){ai,S}=Σi[{1,S}+{ai,∂S/∂ai}]=ΣiΣi[{1,aj}(∂S/∂aj)+{ai,aj}(∂2S/∂ai∂aj)]=ΣiΣi{1,aj}(∂S/∂aj)=-ΣiΣi,k(∂1/∂ak){ak,aj}(∂S/∂aj)=0 となります。1という関数はδ-関数であり,∂1/∂akは汎関数微分ですね。

そこで,kBij=DijはPeq(a)=exp[S()/kB]が解になるための必要十分条件であることがわかりました。

それ故,"揺動力=ゆらぎ"の時間相関関数<Ri(t)Rj(t')>(白色雑音の2Dは揺動力の強さと呼ばれる)と,輸送係数:Lijの間には<Ri(t)Rj(t')>=2kBijδ(t-t')という関係が成り立ちます。

 

これはさらに∫0d(t-t')<Ri(t)Rj(t')>=kBij:すなわち∫0dτ<Ri(t)Rj(t+τ)>=kBijと書き直すことができます。

 

つまり,輸送係数あるいはオンサーガー係数:Lijはゆらぎの時間相関関数で与えられます。この法則を揺動散逸定理と呼びます。

こうして,数式的に表現された形の定理が示されても,これが実際の自然現象において物理的にどのような意味を持つのかを理解しなければ,こうした定理の重要性を認識することはできません。

そこで,1例として熱伝導に適用してみます。

平均流速がゼロの1成分流体中の熱伝導方程式はeを単位質量当たりの内部エネルギーとして∂(ρe)/∂t+∇q=0 で与えられます。

 

ここでqは熱流であり熱拡散の線形近似モデルではq=λ∇(1/T)=-κ∇Tと表わされます。κ=λ/T2は熱伝導率と呼ばれています。

 

この現象論的方程式に対して,これのランジュバン方程式は∂(ρe)/∂t+∇q=-∇(,t)となります。ただし(,t)は熱流qのゆらぎです。

qは熱流ですから,流体の局所流速を(,t)とすると,q=ρe=ρe(d/dt)です。

 

示量変数の1つとして=ρeとすると,の不可逆変化部分に対する表式dai/dt=Σjij(∂S/∂aj)+Ri(t)は,dai/dt=d(ρeri)/dt~(Jq)i=Σjij[∂S/∂(ρerj)] +Ri(t)と書けますが,平衡状態ではρVde=TdS-PdVなので体積一定(dV=0 )ならS=ρeV/T+(定数)です。

それ故,結局(Jq)i~ΣjijV(∂/∂rj)(1/T)となります。これをq=λ∇(1/T):すなわち(Jq)i=λ(∂/∂rj)(1/T)と比較すると,輸送係数=オンサーガー係数についてLij=(λ/V)δijという表式が得られます。

そこで揺動散逸定理によると,λは平衡状態における揺動熱流(,t)の時間相関関数によって表現されることになります。すなわち,揺動熱流(,t)の時間相関関数は,∫0dt<qα(,t)qβ(',0)>eq=kB(λ/V)δαβδ(')と書けます。

そして,熱流q(,t)のゆらぎ(,t)を体積Vの対象領域全体で空間積分した(t)=∫dq(,t)という量を定義して,これを時刻tでの"熱流のゆらぎ"と呼ぶことにすれば,熱伝導に対する揺動散逸定理の表現は∫0dt<qα(t)qβ(0)>eq=kBλδαβと書くことができます。

今は,平均流速がゼロの流体を考えており,平衡状態では<q(,t)>eq=0 なので,熱流のゆらぎ(,t)がエネルギー流そのものになります。そこで対流がない流体においての平衡状態では(t)はエネルギー流を空間全体で積分したもののゆらぎとなります。

この例では,輸送係数=オンサーガー係数の一種である熱伝導率がミクロな流れ(t)の時間相関関数で与えられるということが揺動散逸定理からの重要な帰結と言えます。

19日木曜日夜から風邪気味で,症状自体は軽いのですがあまり筆が進みません。

北原和夫 著「非平衡系の統計力学」(岩波書店)

 

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

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

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

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

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

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

2007年7月19日 (木)

労働貴族死す

  労働貴族の長=宮本顕冶もとうとう死んだか。。。

  とりあえず, 合掌!

 東洋思想では生前に何をしたかによらず死ねば皆仏になるらしい

        

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

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

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

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

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

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

2007年7月18日 (水)

ブラウン運動と伊藤積分(10)

 確率積分の話題の続きです。前回はこれを定義しただけですが今回はこれの性質,特に伊藤の公式について説明します。 

次の不等式は,シュヴァルツ(Schwartz)の不等式の一般化です。

(補題9.1):(国田・渡辺の不等式)

  M,N∈2,cとする。

 

  φ(s,ω),ψ(s,ω)は発展的可測で∀t>0 に対して,∫0tφ(s,ω)2d<M>s<∞,∫0tψ(s,ω)2d<N>s<∞ とする。

 

このとき,∀t>0 に対して,|∫0tφ(s,ω)ψ(s,ω)d<M,N>s|≦(∫0tφ(s,ω)2d<M>s)1/2(∫0tψ(s,ω)2d<N>s)1/2 (a.s)である。 (a.s.はalmost surely="ほとんど確実に"です。)

(証明)<M,N>の定義:<M,N>≡(1/4)(<M+N>-<M-N>)と,その性質:<aM,N>=a<M,N>によって,<M>=<M,M>,<N>=<N,N>であり,∀r∈Rに対して<M+rN>=<M>+2r<M,N>+r2<N>となります。

したがって,0≦∫s1s2d<M+rN>s=∫s1s2d<M>s+2r∫s1s2d<M,N>s+r2s1s2d<N>s a.s が全ての実数rに対して成立するようにできます。

よって,実数rの2次式が常に非負でr2の係数が正なので,その判別式は(判別式)≦0 を満足しなければならないことから,不等式|∫s1s2d<M,N>s|2≦(∫s1s2d<M>s)(∫s1s2d<N>s)が得られます。

それ故,φ(s,ω)≡Σiφi(ω)1(ti,ti+1),ψ(s,ω)≡Σiψi(ω)1(ti,ti+1)0に対して,|∫∫0tφ(s,ω)ψ(s,ω)d<M,N>s|=|Σiφi(ωi(ω)∫titi+1d<M,N>s|≦Σii(ω)||ψi(ω)|(∫0td<M>s) 1/2(∫0td<N>s)1/2≦(Σii|2titi+1d<M>s)1/2ii|2titi+1d<N>s)1/2=(∫0tφ(s,ω)2d<M>s)1/2(∫0tψ(s,ω)2d<N>s)1/2 (a.s)となります。

有界なφ(s,ω),ψ(s,ω)については,再掲(補題:8.3):"02(<M>)内で稠密である"によって,上のような0の元で近似することができます。

 

一般のφ,ψについてはまず有界なもので近似して0の元で近似すると極限で命題が成立します。

 

(証明終わり)

(定理9.2) 

(ⅰ) M,N∈2,c,f(s,ω)∈2(<M>)のとき,Xt=∫0tf(s,ω)dMsとすると,<X,N>t=∫0tf(s,ω)d<M,N>sであり,かつ任意のN∈2,cに対して<X,N>t=∫0tf(s,ω)d<M,N>sを満たすX∈2,cでX0=0 なるものはXt=∫0tf(s,ω)dMsに限る。

 

(ⅱ)I(f)=XT=∫0f(s,ω)dMs2(0,T;<M>)→2,c(0,T)は∀Tに対して等距離写像(isometric mapping)である。:すなわちE[|XT|2]=E[∫0|f(s,ω)|2d<M>s]である。

(証明)(ⅰ)f∈2(<M>)に対してfn0でE[∫0t|fn(s,ω)-f(s,ω)|2d<M>s]→ 0  as n→ ∞なるものを取ります。

 

nt=∫0tn(s,ω)dMsとおくと,再掲(補題8.7):"E[X2]=E[<X>T]=E[∫0Tf(s,ω)2d<M>s]である。"によって,E[<Xn-X>t]=E[∫0t|fn(s,ω)-f(s,ω)|2d<M>s] → 0  as n→ ∞,∀tです。

一方,|∫d<M,N>s|≦|∫d<M>s|1/2|∫d<N>s|1/2ですが,E[|<M,N>t|]=E[|∫0td<M,N>s|],E[|<M>t|]=E[|∫0td<M>s|],E[|<N>t|]=E[|∫0td<N>s|]より,一般にE[|<M,N>t|]≦E[|<M>t|]1/2E[|<N>t|]1/2です。

したがって,E[|<Xn,N>t-<X,N>t|]=E[|<Xn-X,N>t|]≦E[|<Xn-X>t|]1/2E[|<N>t|]1/2 → 0  as n →∞ ,∀tが得られます。

また,<Xn,N>t=∫0tn(s,ω)d<M,N>sがfn(s,ω)∈0を具体的に書き下すことによって(補題8.5)の証明と同じようにして言えるので,結局<X,N>t=∫0tf(s,ω)d<M,N>s a.sであることがわかります。

 後半の,一意性については,<X~,N>t=∫0tf(s,ω)d<M,N>s a.s ∀N∈2,cとすると,<X-X~,N>t=0 a.sですから,特にN=X-X~とすると<X-X~>t=0 a.s,あるいはE[|Xt-X~t|2]=0 よりXt=X~t a.sとなることから示されます。

(ⅱ)(ⅰ)よりXt=∫0tf(s,ω)dMsとおけば,<X>T=<X,X>T=∫0Tf(s,ω)d<M,X>sであり,<M,X>s=∫0sf(u,ω)d<M>uですからd<M,X>s=f(s,ω)d<M>sなので,<X>T=∫0Tf(s,ω) 2d<M>s=∫0T|f(s,ω)|2d<M>sとなります。

 

 これとE[|XT|2]=E[<X>T]を合わせると,E[|XT|2]=E[∫0T|f(s,ω)|2d<M>s]が得られます。(証明終わり)

(系9.3)(ⅰ)M∈2,c,f,g∈2(<M>),a,b∈Rに対して∫0t{af(s,ω)+bg(s,ω)}dMs=a∫0tf(s,ω)dMs+b∫0tg(s,ω)dMs,∀tである。

 

(ⅱ)M,N∈2,c,f∈2(<M>)∩2(<N>),a,b∈Rに対してf∈2(<aM+bN>)で∫0tf(s,ω)d(aM+bN)s=a∫0tf(s,ω)dMs+b∫0tf(s,ω)dNs ,∀tである。

(証明) (ⅰ)N∈2,cのとき,(定理9.2)と二次変分の性質(命題7.8)より,<∫0t(af+bg)dMs,N>=∫0t(af+bg)d<M,N>s=a∫0tfd<M,N>s+b∫0tgd<M,N>s=<a∫0tfdMs,N>+<b∫0tgdMs,N>です。

 

 そしてNは任意なので,例えばN=∫0t(af+bg)dMs-a∫0tfdMs-b∫0tgdMsと取ることによって,∫0t{af(s,ω)+bg(s,ω)}dMs=a∫0tf(s,ω)dMs+b∫0tg(s,ω)dMs ,∀tが得られます。

(ⅱ)f∈2(<aM+bN>)なることについては,再掲(補題9.1):"M,N∈2,cとしφ(s,ω),ψ(s,ω)は発展的可測で,∀t>0 に対して∫0tφ(s,ω)2d<M>s<∞,∫0tψ(s,ω)2d<N>s<∞ とする。このとき∀t>0 に対し|∫0tφ(s,ω)ψ(s,ω)d<M,N>s|≦(∫0tφ(s,ω)2d<M>s)1/2(∫0tψ(s,ω)2d<N>s)1/2 a.sである。"と二次変分の性質からいえます。

 

すなわち,(補題9.1)でφ=ψ=fとおくと,不等式|∫0tf(s,ω)2d<M,N>s|≦(∫0tf(s,ω)2d<M>s)1/2(∫0tf(s,ω)2d<N>s)1/2 a.sが得られます。

 

 これと等式d<aM+bN>s=a2d<M>s+2abd<M,N>s+b2d<N>sから,f∈2(<aM+bN>)なることは自明です。

 そして(定理9.2)と(命題7.8)によれば,∀L∈2,cに対して<∫0tfd(aM+bN)s,L>=∫0tfd<(aM+bN),L>s=a∫0tfd<M,L>s+b∫0tfd<N,L>s=<a∫0tfdM+b∫0tfdN,L>が成立するので命題が成り立つことは明らかです。

 

 (証明終わり)

 次に,局所マルチンゲールに対して確率積分を定義します。 

(定義10.1):(局所マルチンゲールに対する確率積分)

M∈,locとし,Mに対しP(∫0Tf(s,ω)2d<M>s<∞)=1,∀Tを満たす発展的可測過程f(s,ω)について確率積分を定義する。

M∈,locよりσn↑∞ a.sとなる停止時刻の列{σn}が存在して,Mt∧σn2,cが成立する。そしてτn(ω)≡n∧inf{t≧0:∫0tf(s,ω)2d<M>s≧n}とおけばτn↑∞ a.sである。

 

 そこでρn≡τn∧σnとしMnt≡Mt∧ρn,fn(t,ω)≡f(t,ω)1{t≦ρn}とおけば,Mnt2,c,fn(t,ω)∈2(<Mn>)なので確率積分I(fn)がI(fn)(t,ω)=∫0tn(s,ω)dMnsと定義できる。

 

 これはI(fn)(t,ω)=I(fm)(t,ω),0≦t≦ρn,n≦mとなることがわかる。

 そして,n→ ∞ に対してσn → ∞ なので,I(f)(t,ω)≡I(fn)(t,ω), 0≦t≦ρnとすれば,n→ ∞ の極限で∀tについてI(f)が定義できる。

  

 このI(f)をfの確率積分という。

(定理9.2)を局所マルチンゲールに対して書き直したものはM∈c,loc ,f(s,ω)をP(∫0Tf(s,ω)2d<M>s<∞)=1,∀Tを満たす発展的可測過程とする。

 

 このときXt=∫0tf(s,ω)dMs,M∈c,locはX0=0 で<X,N>t=∫0tf(s,ω)d<M,N>s∀t a.sを満たす唯1の元である。

 

 となります。

 

 これの証明は,tをt∧ρnにおきかえると,(定理9.2)の証明と同じなので省略します。

次に,確率積分が普通の積分と同様な概念を表現するものであるということを示す重要な性質を証明します。

 "f(s,ω)がtに適合して左連続であるとする。このとき,分割Δ:0=s0<s1<..<snをとれば,∫0tf(s,ω)dMs=P-lim|Δ|→0Σif(si)(Msi+1-Msi)となる。ただし,P-lim は確率収束極限の意味である。という命題を証明します。

(証明)fが有界でf∈2,cのとき,分割Δ:0=t0<t1<..<tnに対してfΔ(s,ω)≡Σif(ti,ω)1(ti,ti+1)(s),fΔ(0,ω)≡f(0,ω)とおくと,左連続性によってfΔ(s,ω)→ f(s,ω) as Δ→ 0,∀s,ωです。

Δt=∫0tΔ(s,ω)dMs=Σif(ti,ω)(Mti+1-Mti)であり,E[∫0t|fΔ(s,ω)-f(s,ω)|2d<M>s]→ 0 よりE[|XΔt-Xt|2] → 0 as Δ→ 0,∀tです。

一般の局所マルチンゲールの場合は(定義10.1)のτnを使って局所化すればいいだけです。

 

(証明終わり)

ここで微積分学における合成関数の微分法則(連鎖公式)の確率解析版とされる伊藤の公式を示すことにします。

 

その応用は極めて広いものです。

(定義11.1)確率過程X={Xt}がXt=X0+Mt+At,M∈loc,A∈(増加過程:に属する元の差で表わされる過程)で,X00-可測関数,ただし,M0=0 a.s,A0=0 a.sと表わされるとき,{Xt}は半マルチンゲールであるという。

また,RN値確率過程が半マルチンゲールであるとは,各成分が半マルチンゲールのときをいう。

 

さらに,t,tが連続確率過程であれば{t}は連続な半マルチンゲールであるという。

(定理11.2)(伊藤の公式)

t(X1t,X2t,..,XNt)を連続な半マルチンゲールとする。すなわち,Xit=Xi0+Mit+Ait (i=1,2,..,N)とする。

 

f∈2(RN)のとき,f(t)は連続な半マルチンゲールであり,f(t)-f(0)=Σi=1N0tif(s)dMis+Σi=1N0tif(s)dAsi+(1/2)Σi,j=1N0tijf(s)d<Mi,Mjsと書くことができる。

 

ただし,Dif≡∂f/∂xi,Dijf≡∂2f/∂xi∂xj(i,j=1,2,..N)である。

(証明)τn=inf{t≧0:|0|>n,or|t|>n,|t|>n}({ }≠φのとき),τn=∞ ({ }=φのとき);とおきます。このとき,n→ ∞ に対してτn→ ∞ a.sなので,Xt∧τn について,この等式を証明すれば十分です。

したがって,0,t,tは全て有界,f,Dif,Dijfも全て有界,かつ一様連続であるとしてかまいません。それ故,ある定数Kがあって|f|+Σi|Dif|+Σi,j |Dijf|≦Kと書くことができます。

Δ:0=t0<t1<..<tn=tを [0,t]の分割とします。

このとき,テイラー展開の定理により,f(t)-f(0)=Σk=1n-1(f(tk+1)-f(tk))=Σk=1n-1Σi=1Nif(tk)(Xitk+1-Xitk)+(1/2)Σk=1n-1Σi,j=1Nijf(ξk)(Xitk+1-Xitk)(Xjtk+1-Xjtk),ξktk+θk(tk+1tk), 0≦θk≦1と表現できます。

右辺の第1項は|Δ|→0 のとき,確率積分の定義により,Σi=1N0tif(s)dMis s+Σi=1N0tif(s)dAsiに確率収束します。

一方,右辺の第2項は,(1/2)Σk=1n-1Σi,j=1Nijf(ξk)(Mitk+1-Mitk)(Mjtk+1-Mjtk)+Σk=1n-1Σi,j=1Nijf(ξk)(Mitk+1-Mitk)(Ajtk+1-Ajtk)+(1/2)Σk=1n-1Σi,j=1Nijf(ξk)(Aitk+1-Aitk)(Ajtk+1-Ajtk)≡I1+I2+I3となります。

 

より,,,+,cとすると,|Σk=1n-1(tk+1tk)|≦|t|+|t|で,sup|tk+1tk|→ 0,sup|tk+1tk|→ 0 なので,|I2|≦Ksup k|tk+1tk|(|t|+|t|)→ 0 as |Δ|→ 0 a.s,かつ|I3|≦Ksup|tk+1tk|(|t|+|t|→ 0 as |Δ|→ 0 a.sです。

よって,後はI1(1/2)Σi,j=1N0tijf(s)d<Mi,Mj> as |Δ|→ 0 a.sを示すことができればいいことになります。

 

通常の積分と微分との関係では2次の微小量を総和して積分しても1次の微小量になるだけで,|Δ|→ 0 では消えてしまうのでこうした2階導関数の項は出現しません。

実際,上に示したように有界変動の関数と有界な単調増加関数の差で表わせるので,その2次の量を積分するとき,その寄与はゼロになります。

 

しかし,がブラウン運動のような過程である場合には有界変動の関数ではなくて,その長さが確率1で無限大になることは,ずいぶん前の記事で述べましたが,その場合のの2次の量を積分するとき,その寄与はゼロではなくて有限です。

1'=(1/2)Σk=1n-1Σi,j=1Nijf(tk)(Mitk+1-Mitk)(Mjtk+1-Mjtk)とおくと,|I1-I1'|≦(1/2)Σk=1n-1Σi,j=1N|Dijf(ξk)-Dijf(tk)||Mitk+1-Mitk||Mjtk+1-Mjtk|≦(1/2)supk|Dijf(ξk)-Dijf(tk)|[Σn,mt(Mn,Δ)1/2t(Mm,Δ)1/2]です。

よって,E[|I1-I1'|]≦(1/2)Σn,m=1NE[supk|Dijf(ξk)-Dijf(tk)|Qt(Mn,Δ)1/2t(Mm,Δ)1/2]≦(1/2)Σn,m=1NE[supk|Dijf(ξk)-Dijf(tk)|2]1/2(E[Qt(Mn,Δ)Qt(Mm,Δ)])1/2≦(1/2)E[supk|Dijf(ξk)-Dijf(tk)|2]1/2n,m=1NE[Qt(Mn,Δ)2]1/4E[Qt(Mm,Δ)2]1/4)です。

 

再掲(補題7.3):"M={Mt}∈4,cとするとき分割Δに依らない定数cが存在してE[Qt(M;Δ)2]≦cが成り立つ。"とE[supk|Dijf(ξk)-Dijf(tk)|2]1/2→ 0 as |Δ|→ 0 によって,右辺→ 0 as |Δ|→ 0 が得られます。

また,I1"=(1/2)Σk=1n-1Σi,j=1Nijf(tk)(<Mi.Mjtk+1-<Mi.Mjtk)とおくとE[|I1'-I1"|2]=(1/4)Σk=1 n-1E[|Σi,j=1Nijf(tk){(Mitk+1-Mitk)(Mjtk+1-Mjtk)-∫tktk+1d<Mi.Mjs|2]≦Σk=1 n-1(K2/4)E[|tk+1tk|4+Σi,j=1N(<Mi.Mjtk+1-<Mi.Mjtk)2]≦(K2/4)E[supk|tk+1tk|2Σi=1Nt(Mi,Δ)]+(K2/4)Σi,j=1NE[supk|<Mi.Mjtk+1-<Mi.Mjtk|(<Mi.Mjtk+1-<Mi.Mjtk)] → 0 as |Δ|→ 0 です。

 

すなわち,E[|I1'-I1"|2]→ 0 as |Δ|→ 0 が得られます。

一方,1"=(1/2)Σk=1n-1Σi,j=1Ntktk+1ijf(tk)d<Mi.Mjs→(1/2)Σi,j=1N0tijf(s)d<Mi.Mjs as |Δ|→ 0 ですから,結局,I1→(1/2)Σi,j=1N0tijf(s)d<Mi.Mjs in L(Ω) as |Δ|→ 0 です。

以上から,∀tに対して,ほとんどいたるところで,f(t)-f(0)=Σi=1N0tif(s)dMis+Σi=1N0tif(s)dAsi+(1/2)Σi,j=1N0tijf(s)d<Mi,Mjsが成立し,両辺はtについて連続なので,この等式は確率1で∀tに対して成立します。

 

(証明終わり)

例として原点から出発する1次元ブラウン運動をXt=Btとし,f(x)=xnとすると,伊藤の公式から,Btn=n∫0tsn-1dBs+{n(n-1)/2}∫0tsn-2dsです。

 

特にn=2とすれば,Bt2=2∫0tsdBs+tです。

 

このことはBt2-tがマルチンゲールであることを示しています。これは,マルチンゲールの確率積分が,その定義によって常にマルチンゲールであるからです。

こうして,1次元ブラウン運動の確率積分が通常の関数についての積分公式t22∫0txdxとは食い違うことが明確にわかります。

 次の例として,At=-t/2としXt=Bt+At=,f(x)=exp(x)とするとき,伊藤の公式から,exp(Xt)-exp(X0)=exp(Bt-t/2)-exp(B0)=∫0t exp(Bs-s/2)dBs+∫0t exp(Bs-s/2)d(-s/2)+(1/2)∫0t exp(Bs-s/2)dsです。

 

したがって,exp(Bt-t/2)=1+∫0t exp(Bs-s/2)dBsです。

tはマルチンゲールなので,もちろんその確率積分exp(Bt-t/2)もマルチンゲールです。

 

また,E[exp{iξ(Bt-Bs)}]=exp{-(t-s)ξ2/2}より,s=0 ,ξ=-2iとおいて,E[exp(2Bt)]=exp(2t),すなわち,E[exp(2Bt-t)]=exp(t)です。

 

それ故,exp(Bt-t/2)は2乗可積分マルチンゲールです。

さらにt(B1t,B2t,..,BNt)をN次元ブラウン運動とし,f()={(x1)2(2)2+..+(xN)2}m/2とするとき,m≧2ならf(t)-f(0)=mΣi=1N0tis|s|m-2dBis+(1/2)∫0t{Nm+m(m-2)}|s|m-2dsです。

 

また,N≧3でm=2-N,σn=inf{t:|t|≦1/n}とするとき,f(t∧σn)-f(0)=(2-N)Σi=1N0t∧σnis|s|-NdBisですから,tの出発点が 0 でないとき:P(t0)=1,00 なるときも含めて,P(σn↑∞ as n→∞)=1 なので,f(t)-f(0)は局所マルチンゲールです。

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

参考文献:長井英生 著「確率微分方程式」(共立出版)

 

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

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

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

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

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

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

目をそらす姑息な手段

 しかし,この赤城という農相の姑息な手段にもあきれますね。どんな参謀がついているのでしょうか?まあ,それに踊らされる方もどうかと思いますね。

 顔に絆創膏というなんともお粗末なカムフラージュ,本題の方の追求からまわりの目をそらすために,あえて突っ込まれても痛くない別の目立った疑惑を無理矢理作り上げる,という昔からよくある手法ですね。

 「不利なときは戦線を拡大せよ。」という将棋の戦法にも通じますが,これまでも何か不利なことがあってもすぐに別の問題があって興味がそらされてしまうというのがこのところの政界では日常茶飯事ですね。

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

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

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

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

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

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

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年7月12日 (木)

ブラウン運動と伊藤積分(9)

 次のステップとして確率積分という主題に移ります。 

以下では,フィルター付確率空間(Ω,,P;t)はtが右連続であり,かつ次の条件を満たしているものとします。すなわち,≡{A∈Ω|∃B∈:A⊂B,P(B)=0}⊂0 であるとします。

2(0,T)≡{M:E[|MT|2]<∞;{Mt}t∈[0,T]は右連続マルチンゲール},2,c(0,T)≡{M∈2(0,T):{Mt}t∈[0,T]は連続}と置きます。

 

また,特にσ-加法族の増大系tを強調する場合には,2(0,T;t),2(t),2,c(0,T;t),2,c(t)のように表わします。

(系7.6)によれば,∀M∈2,cに対し,ある<M>∈+,cがあって,Mt2-<M>tはマルチンゲールになっています。