ラプラス方程式→連立方程式

質問

楕円型ラプラス方程式を連立方程式に近似して解きたいと思っています。
そこで、いろいろ調べてみたのですが、どうしても、
連立方程式の基本形　Ax=b　のうち、ラプラス方程式の差分近似方程式のどの部分がどの文字に当たるのかが理解できません。

例として以下のラプラス方程式を連立方程式Ax=bの形そのものに、直していただけないでしょうか。

(∂^2 * u) / (∂^2 * x)　+　(∂^2 * u) / (∂^2 * y) = 0 (0<x<1 , 0<y<1)
u(x,0)=sinπx , u(x,1)=0　(0≦x≦1)
u(0,y)=u(1,y)　(0≦y≦1)

近似方程式は
u(x_i,y_j) =
(u(x_i+h , y_j) + u(x_i-h , y_j) + u(x_i , y_j+h) + u(x_i , y_j-h)) / 4

です。

よろしくお願いします。

回答No.1

stomachman · ベストアンサー

hってのは、
x_i = hi, y_j = hj, i∈{0,1,...,N}, j∈{0,1,...,N}, N=1/h
という意味なのかな。ま、そうだとしましょう。

　連立方程式の変数をとりあえずu[i,j]と書くことにして、
u[i,j] = u(x_i, y_j)　　(i∈{0,1,...,N-1}, j∈{0,1,...,N})
としますと、u[i,j]を１列に並べて作った(N(N+1))次元ベクトルがxです。並べ方はどんな風でも良い。デタラメに並べたって構いません。ですが、
(1)　I(k) = (x[k]に対応するu[i,j]のi)
(2)　J(k) = (x[k]に対応するu[i,j]のj)
(3)　K(i,j) = (u[i,j]に対応するx[k]のk)
を表にするなり、計算アルゴリズムを作るなりしておかないと、分かんなくなります。
　なお、u[N,j] が現れていません。というのは、u[N,j]= u[0,j]なので、変数にしなかったのです。（もちろん、u[N,j]をu[0,j]とは別の変数とみなした上で、方程式u[N,j]= u[0,j] (j∈{0,1,...,N})を追加しても構わないのですが。）

　A, bに要素がどういう順番で並ぶかは、u[i,j]を１列に並べてxを作る際に、どう並べたかによって決まります。方程式は全部一次式で
u[i,0]=sin(πhi) (i∈{0,1,...,N-1})
u[i,N]=0 (i∈{0,1,...,N-1})
の(2N)本と、
u[i,j]-(1/4)u[i+1,j]-(1/4)u[i-1,j]-(1/4)u[i,j+1]-(1/4)u[i,j-1]=0 　(i∈{1,...,N-2}, j∈{1,...,N-1})
の((N-2)(N-1))本、それに、iの両端（i=0とi=N-1）においてu[0,j]=u[N,j]を代入した
u[0,j]-(1/4)u[1,j]-(1/4)u[N-1,j]-(1/4)u[0,j+1]-(1/4)u[0,j-1]=0 　( j∈{1,...,N-1})
u[N-1,j]-(1/4)u[0,j]-(1/4)u[N-2,j]-(1/4)u[N-1,j+1]-(1/4)u[N-1,j-1]=0 　(j∈{1,...,N-1})
の(2(N-1))本、合計で(N(N+1))本あるので変数の個数と同じ。

　従って、係数の行列Aは(N(N+1))行(N(N+1))列の正方行列になりますね。対角要素は1で、あと-(1/4)がパラパラ入っている以外は、どの要素も0です。
　定数のベクトルbは(N(N+1))次元で、これまたほとんどの要素が0だけど、u[i,0]=sin(πhi) に当たる部分だけ0でない値になっています。

お礼

理解できました。
ありがとうございました！