« 次元控除法によるポアソン方程式の直接数値解法(1) | トップページ | 次元控除法によるポアソン方程式の直接数値解法(補遺) »

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/「健康商品の店 タクザイ」

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

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

|

« 次元控除法によるポアソン方程式の直接数値解法(1) | トップページ | 次元控除法によるポアソン方程式の直接数値解法(補遺) »

308. 微分方程式」カテゴリの記事

311 .数値計算・調和解析・離散数学」カテゴリの記事

コメント

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: 次元控除法によるポアソン方程式の直接数値解法(2):

« 次元控除法によるポアソン方程式の直接数値解法(1) | トップページ | 次元控除法によるポアソン方程式の直接数値解法(補遺) »