« 超幾何微分方程式の解と接続公式 | トップページ | 次元控除法によるポアソン方程式の直接数値解法(2) »

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) 他

|

« 超幾何微分方程式の解と接続公式 | トップページ | 次元控除法によるポアソン方程式の直接数値解法(2) »

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

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

コメント

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

トラックバック


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

« 超幾何微分方程式の解と接続公式 | トップページ | 次元控除法によるポアソン方程式の直接数値解法(2) »