1 楕円型方程式

楕円型方程式は、典型的にはラプラス方程式

$$\nabla^2 u = 0$$ (1)

またはポアソン方程式

$$\nabla^2 u = -4\pi G \rho$$ (2)

のようなものである。上のポアソン方程式はニュートン重力の場合で、 $G$ は重力定数、$\rho$ は物質の質量密度（単位体積あたりの質量）である。 ここで、 $\rho$ はとりあえず与えられた空間座標だけの関数とする。空間1次元ではラプラシアンがただの2階微分演算子になって しまい、解くべき方程式は常微分方程式ということになる。

実際にポアソン方程式（とさらに状態方程式を連立させてでてくるレーン・エ ムデン方程式）を球対称の場合に解くというのは I のほうでやったと思うの で、ここでは空間2次元の場合を考える。3次元も考え方は同様である。

工学的な応用では、楕円型方程式の応用として重要なのはラプラス方程式の境 界値問題である。これは、構造解析、電場の解析などあらゆるところで現れる。 特に複雑な形状の機械部品とか機械全体の中での応力場の解析といったものが 重要な応用ということになる。

が、天文・物理のいろんな問題では、複雑な境界値問題というのはあんまり大事では ない。天体は大抵重力でまとまっているわけで、固定された表面（境界）とい うものが与えられているとは限らないからである。もちろん、数値計算の方法 によってはそういうものが出てくることもあるが、とりあえず今日は比較的に 単純な境界条件だけをかんがえることにしよう。

ラプラス方程式やポアソン方程式は、基本的には拡散方程式から時間微分の項 を落したものと考えることができる。つまり、なんらかの定常状態を記述して いるわけである。多くの応用で、実際にそのような定常状態を求めているわけ であるが、数値計算という観点からすると、ラプラス方程式やポアソン方程式 に適用できる計算法は実は拡散方程式や波動方程式に陰解法を適用する時にも 使える。

これは、前回見たように陰解法では空間方向の差分方程式を解く必要がでてく るからである。解くべき方程式は、ラプラス方程式やポアソン方程式の差分化 からでてくるものと基本的には同じである。

というわけで、今日は、単純な2次元ポアソン方程式の境界値問題を例に、計 算方法について考えていく。

2次元なので、ラプラシアンは

$$\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$$ (3)

である。例によって変数の方を $u_{ij}$ と書き、これが $u(x_i,y_j)$ の近 似解としよう。また、前回と同じ2次精度の中心差分を使うとすれば、ラプラ シアンの差分近似は以下で与えられる。

$$\frac{u_{i,j+1}+u_{i,j-1}+u_{i+1,j}+u_{i-1,j}-4u_{ij}} {\Delta x^2}$$ (4)

但し、ここでは $\Delta x = \Delta y$ とした。いま、解くべき問題を長方 形領域でのポアソン方程式の境界値問題

$$\nabla^2 u(x,y) = \rho(x,y), \quad(x_0\le x \le x_1, y_0 \le y< \le y_1)$$ (5)

$$u(x_0,y) = u(x_1,y) = u(x,y_0) = u(x,y_1) = 0$$ (6)

ということにすれば（以下、面倒なので $-4\pi G$ の項は 1 とする）、 これは2次元配列u[nx][ny]についての線形連立方程式ということになる。 ここで

$$n_x = (x_1-x_0)/\Delta x + 1$$ (7)

$$n_y = (y_1-y_0)/\Delta x + 1$$ (8)

であり、境界上では $u_{ij}=0$, それ 以外では

$$\frac{u_{i,j+1}+u_{i,j-1}+u_{i+1,j}+u_{i-1,j}-4u_{ij}} {\Delta x^2} = \rho_{ij}$$ (9)

となる。