次元控除法によるポアソン方程式の直接数値解法(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);
z0=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,k,φNx+1,j,k,φi,0,k,φi,Ny+1,k,φi,j,0,φi,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)となります。
Φk≡t(φ1,k,φ2,k,..,φNy,k),
φj,k≡t(φ1,j,k,φ2,j,k,..,φNx,j,k),
Ρ'k≡t(ρ'1,k,ρ'2,k,..,ρ'Ny,k)
と定義します。
ρ'j,k≡t(ρ1,j,k+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/hz,ρ2,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)によって,ΛNx=tTNxATNx (ΛNxは対角成分
がλ1,λ2,..,λ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)
この得られた対角行列をΩとし,Pk≡tQ~Φk,Rk≡tQ~Ρ'k
(k=1,2,..,Nz)とすれば,(3-4)のPoisson方程式は,
Pk+1-2Pk+Pk-1+ΩPk=Rk(k=2,3,..,Nz-1) (3-18a),
P2-P1+ΩP1=R1 (3-18b),
-PNz+PNz-1+ΩPNz=RNz (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次元の連立1次方程式系であり,係数行列は,i=j=1のとき非
正則,その他の場合は正則ですから,2次元の場合と同様,容易に
解けます。
1次元の連立1次方程式を解いてPkが得られれば,
それに対して,Φk=Q~Pkとすれば,最終的な解Φ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(TN)lm=(1/N)1/2 (m=1),
(TN)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としすると,
既に示したように,固有値問題Ax=λx(x≠0)の解として,
N個の異なる固有値λ=λm=-4sin2{(m-1)π/(2N)}
(m=1,2,..,N)が得られます。
そして,固有値λに属する固有ベクトルxは,
θ=θm=(m-1)π/N (λ=λm=-4sin2(θm/2) )として,
x=t(x1,x2,..,xN)の成分xkが,
xk=a[exp(ikθ)+exp{-i(k-1)θ}}
=2aexp(-iθ/2)cos{(2k-1)θ/2}なる形で与えられること
も,既に前記事の(注2)で見た通りです。
そこで,λmに属する固有ベクトルを改めてxm=t(x1m,x2m,..,xNm)
とおけば,xlm=cmcos{(2l-1)(m-1)π/(2N)}です。
そして,固有値λ1,λ2,..,λN はすべて異なるので,
m≠nなら実内積(xmxn)は直交します,つまり(xmxn)=txmxn=0
を満たします。
一方,(xmxm)=txmxm=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)とおけば,(xmxn)=txmxn=δmnと正規直交化
されます。
そこで,N×N行列TNで成分が(TN)lm=xlmとなるものを与えると,
このTN は1つの直交行列で,Aは直交変換tTNATN=Λ(Λは対角
成分がλ1,λ2,..,λNの対角行列)によって対角化されます。
(以上)
参考文献:1)Masako Ogura:”A direct solution of Poisson's equation by dimention reduction method.”Journal of the Meteorogical Soc. Japan. Vol.47. No.4(1969).pp319-322(日本気象学会誌、第47巻4号,1969) 2)矢嶋信夫、野木達男 著「発展方程式の数値解法」(岩波書店) (1977)他
http://folomy.jp/heart/「folomy 物理フォーラム」サブマネージャ ー
人気blogランキングへ ← クリックして投票してください。(1クリック=1投票です。1人1日1投票しかできません。)
http://homepage2.nifty.com/toshis-kaiga-auction/「健康商品の店 タクザイ」
← クリックして投票してください。(ブログ村科学ブログランキング)
![]() 物理学 |
| 固定リンク
「308. 微分方程式」カテゴリの記事
- 積分方程式(2)(2009.09.11)
- 積分方程式(1)(導入)(2009.08.30)
- 水の波(8)(有限振幅の波:非線形波3)(2009.07.24)
- 水の波(7)(有限振幅の波:非線形波2)(2009.07.17)
- 水の波(6)(有限振幅の波:非線形波1)(2009.07.15)
「311 .数値計算・調和解析・離散数学」カテゴリの記事
- 準線形計画法(ネットワーク輸送問題)(3)(2007.07.16)
- 準線形計画法(ネットワーク輸送問題)(2)(2007.07.15)
- 準線形計画法(ネットワーク輸送問題)(1)(2007.07.13)
- 次元控除法によるポアソン方程式の直接数値解法(補遺)(2007.05.21)
- 次元控除法によるポアソン方程式の直接数値解法(2)(2007.05.20)
この記事へのコメントは終了しました。
コメント