楕円型偏微分方程式

ポアソン方程式
φ[i+1,j]-2φ[i,j]+φ[i-1,j]/(Δx)^2
+φ[i,j+1]-2φ[i,j]+φ[i,j-1]/(Δy)^2
=-13π^2sin3πxsin2πy
の解の求め方が分かりません。

解法の手順のご教授お願いします。
最終的には連立一次方程式の形に持って行き、
ガウスの消去法で解けるようにしたいです。
そこまで辿り着くのにどのように解いていけばいいのか分からないです…。
まずは何をすればいいのか？
分点の数を決めればいいのか？格子点の数を決めればいいのか？

No.2 ベストアンサー 回答者： inara 回答日時： 2007/04/14 14:49

差分法による定常2次元ポアソン方程式の解法は、ほとんどの数値計算のテキストに出ているので、Webにもいろいろありますが[1]、エッセンスだけまとめておきます。

参考資料にはプログラムも載っているので、手っ取り早く計算したいだけならそれを使ってください。

フローは以下のとおりです。
【１】　支配方程式と境界条件を差分方程式に変換する
【２】　連立方程式を作る
【３】　連立方程式を解く

【1】　差分方程式の導出
一例として簡単な境界条件（周辺値＝一定：ディリクレ境界条件）での解き方を紹介します。
関数 φ = φ( x , y ) は
支配方程式：∂^2(φ)/∂x^2 + ∂^2(φ)/∂y^2 =- f ( x,y ) ; 0 ≦ x ≦ 1, 0 ≦ y ≦ 1 --- (1)
境界条件：境界上で φ =a　--- (2)
を満足するとします。質問の問題は f(x,y) = 13*π^2*sin(3*π*x)*sin(2*π*y) ですが、一般的に f(x,y)と書いておきます。

２階偏導関数の差分表式は５点近似と９点近似があります[4]が、簡単な５点近似を使って式 [1] を差分方程式に変換すると
{ φ[ i+1 , j ] - 2*φ[i,j] + φ[ i-1 , j ] } / (Δx)^2 + { φ[ i , j+1 ] -2*φ[i,j] + φ[ i , j-1] } / (Δy)^2 = f(x[i] , y[i] ) --- (3)
となります（ここまではVTR250さんの質問に書かれています）。関数 f( x , y ) の x 依存性と y 依存性が著しく異なっていない場合（ f( x , y ) = 13*π^2*sin(3*π*x)*sin(2*π*y) は、x方向もy方向も似た形状なのでOK ）、Δx = Δy = h とおけば式 (3) は
φ[ i+1 , j ] + φ[ i , j+1 ] + φ[ i-1 , j ] + φ[ i , j-1 ] -4*φ[ i , j ] = -h^2*f( i*h, j*h ) ( i = 1, 2, ... , N-1 ; j = 1, 2, ... , N-1 ) --- (4)
となります。一方、境界条件 [2] は以下のようになります。
φ[ i , 0 ] = φ[ i , N ] = a 　( i = 0, 1, 2, ... , N ) 、φ[ 0 , j ] = φ[ N , j ] = a 　( j = 0, 1, 2, ... , N ) --- (5)

【2】　連立方程式を作る
【1】の結果をまとめます
φ[ i+1 , j ] + φ[ i , j+1 ] + φ[ i-1 , j ] + φ[ i , j-1 ] -4*φ[ i , j ] = -h^2*f( i*h, j*h ) ( i = 1, 2, ... , N-1 ; j = 1, 2, ... , N-1 )
φ[ i , 0 ] = φ[ i , N ] = a 　 （ i = 0, 1, 2, ... , N )
φ[ 0 , j ] = φ[ N , j ] = a 　( j = 0, 1, 2, ... , N )

これをベクトル表示 A*φ = B の形（連立方程式）にすればφ が求められます。Aは係数行列（ (N-1)行(N-1)列 ）で、資料[3]の式(4.22）のような形になります。φ は解ベクトル（列ベクトル）で
φ = ( φ[1,1], φ[2,1], ... , φ[1,2], φ[2,2], ... , φ[1, i ], φ[2, i ], ... , φ[1, N-1 ], φ[2, N-1], ... , φ[N-1, N-1 ] )
という形になります。解ベクトルのφ の[ ] 内部は 1 から N-1 までですが、[ ] 内部が0やNになるのは境界上なので、この場合解くまでもなく、最初から値が決まっている（φ [ ] = a )ため除外しています。境界条件がノイマン条件の場合、境界内部のφ の値によって、境界上の値が変わってきます（資料[3]はノイマン型境界条件のポアソン方程式の差分解法です）。解ベクトルの並べ方は任意ですが、これはナチュラルオーダリングと呼ばれる並べ方です。Bは項ベクトル（列ベクトル）で、この場合
B = ( -h^2*f( h, h ) , -h^2*f( 2*h, h ), ... , -h^2*f( h, 2*h ), -h^2*f( 2*h, 2*h ), ... , -h^2*f[1, i ],-h^2*f[2, i ], ... , -h^2*f[1, N-1 ], -h^2*f[2, N-1], ... , -h^2*f[N-1, N-1 ] )
となります。

【3】　連立方程式を解く
【2】の結果をまとめると
A*φ = B
φ = ( φ[1,1], φ[2,1], ... , φ[1,2], φ[2,2], ... , φ[1, i ], φ[2, i ], ... , φ[1, N-1 ], φ[2, N-1], ... , φ[N-1, N-1 ] )
A = (　） --- 資料[3]の式(4.22）の形
B = ( -h^2*f( h, h ) , -h^2*f( 2*h, h ), ... , -h^2*f( h, 2*h ), -h^2*f( 2*h, 2*h ), ... , -h^2*f[1, i ],-h^2*f[2, i ], ... , -h^2*f[1, N-1 ], -h^2*f[2, N-1], ... , -h^2*f[N-1, N-1 ] )
となります。これをφについて解く方法は、資料[2]に出ていますのでここで改めて書きませんが、計算を進めていくと、行列Aが上三角行列になるので、解ベクトルの最下位φ[N-1, N-1 ] から順番に値が求まっていくというものです（資料[3]の例で言えば、最初にx3が求まり、その値を使ってx2が求められる）。φ が求められたら、解ベクトルの並べ方から逆算して、１次元の配列 φ [ 1～(N-1)^2 ] を２次元配列 φ [x,y] に変換します。なお、Gaussの消去法は対角成分(i=j)がゼロになると計算できないので、他の解法も紹介しておきます[5]。

【参考資料】
[1] 数値計算法演習講義ノート（３３ページ：２次元ポアソン方程式）　 http://www.kitasato-u.ac.jp/sci/resea/buturi/his …
[1] CG法によるポアソン方程式の差分解法 http://www.na.cse.nagoya-u.ac.jp/~reiji/lect/alg …
[2] Gaussの消去法（アルゴリズム・Fortranコード）　http://www.na.cse.nagoya-u.ac.jp/~yamamoto/lectu …　の３ページ、または http://polaris.s.kanazawa-u.ac.jp/chuou/chuou.pdf　の８ページ
[3] ２次元 ポアソン方程式の係数行列の構造（p.35-38) 　http://www.coins.tsukuba.ac.jp/~itosho/NC05/chap …
[4] ２階偏導関数の５点近似と９点近似（p.23-24の式8.2と8.3）　http://toshi1.civil.saga-u.ac.jp/ohgushik/sabun1 …
[5] 改良された連立方程式の解法　http://www.na.cse.nagoya-u.ac.jp/~yamamoto/lectu …

この回答へのお礼

すごく丁寧な解説、大変ありがとうございます！
フローチャートだけではなく、参考資料まで…。
感謝感激です！

お礼日時：2007/04/14 18:51

No.1 回答者： komimasaH 回答日時： 2007/04/14 02:13

解があるんですか

この回答への補足

すいません。自分、本当に頭悪くて…。
解を求めるかどうか分かりませんが、
φ[i,j]：未知数＝格子点の数を連立一次方程式に持って行き、ガウスの消去法を用いて導け…と言われたのですが…。

補足日時：2007/04/14 13:03