next up previous contents index
Next: 放物線型偏微分方程式 Up: 偏微分方程式 Previous: 3つの偏微分方程式   目次   索引

楕円型偏微分方程式

楕円型偏微分方程式として,Poisson方程式の境界値問題

$\displaystyle \frac{\partial^{2}u}{\partial x^{2}}(x,y) + \frac{\partial^{2}u}{\partial y^{2}}(x,y) = f(x,y)$ (50)

を平面 $ R = \{(x,y) : a < x < b, c < y < d\}$で考える。ただし, $ u(x,y) = g(x,y), \ (x,y) \in \partial R$. もし,$ f$$ g$$ R$上で連続ならば,この方程式は一意の解を持つ。

ここで,用いる方法は境界値問題への差分法の応用である。つまり,微分方程式における微分を適当な差分係数(difference quotient)で置き換えて解く方法である。 最初のステップは,整数$ m$$ n$を選び,刻み幅 $ h = (b-a)/n$ $ k = (d-c)/m$を決めることである。次に,区間$ [a,b]$を幅$ h$$ n$等分し,区間$ [c,d]$を幅$ k$$ m$等分する。

図 5.1:

\includegraphics[width=8cm]{NUMERICAL/figure10-1.eps}

$\displaystyle x_{i} = a + ih, \ i = 0,1,2,\ldots,n$

$\displaystyle y_{i} = c + jk, \ j = 0,1,2,\ldots,m$

で与えられる点 $ (x_{i},y_{i})$を通る水平線と垂直線を長方形$ R$に引く。

直線$ x = x_{i}$$ y = y_{j}$格子線(grid lines)で,格子線の交点は格子のメッシュ点(mesh points)である。格子の内部のメッシュ点に対しては,$ x_{i}$の周りの$ x$におけるTaylor級数を用いた中心差分公式

$\displaystyle \frac{\partial^{2}u}{\partial x^{2}}(x_{i},y_{i}) = \frac{u(x_{i+...
...)}{h^{2}} - \frac{h^{2}}{12}\frac{\partial^{4}u}{\partial x^{4}}(\xi_{i},y_{j})$ (51)

を用いる。ただし, $ \xi_{i} \in (x_{i-1},x_{i+1})$. 同様にして,$ y_{i}$の周りの$ y$におけるTaylor級数を用いた中心差分公式

$\displaystyle \frac{\partial^{2}u}{\partial y^{2}}(x_{i},y_{i}) = \frac{u(x_{i+...
...}{k^{2}} - \frac{k^{2}}{12}\frac{\partial^{4}u}{\partial y^{4}}(x_{i},\eta_{j})$ (52)

を用いる。ただし, $ \eta_{j} \in (y_{i-1},y_{i+1})$.

これらの2つの公式を用いて式(5.4)を表すと, $ i = 1,2,\ldots,n-1$ $ j = 1,2,\ldots,m-1$に対して,

    $\displaystyle \frac{u(x_{i+1},y_{j}) - 2u(x_{i},y_{j}) + u(x_{i-1},y_{j})}{h^{2}} + \frac{u(x_{i},y_{j+1}) - 2u(x_{i},y_{j}) + u(x_{i},y_{j-1})}{k^{2}}$  
  $\displaystyle =$ $\displaystyle f(x_{i},y_{j}) + \frac{h^{2}}{12}\frac{\partial^{4} u}{\partial x...
...},y_{j}) + \frac{k^{2}}{12}\frac{\partial^{4}u}{\partial y^{4}}(x_{i},\eta_{j})$  

となる。また,境界条件は
$\displaystyle u(x_{0},y_{j})$ $\displaystyle =$ $\displaystyle g(x_{0},y_{j}), \ u(x_{n},y_{j}) = g(x_{n},y_{j}), \ j = 0,1,\ldots,m$  
$\displaystyle u(x_{i},y_{0})$ $\displaystyle =$ $\displaystyle g(x_{i},y_{0}), \ u(x_{i},y_{m}) = g(x_{i},y_{m}), \ i = 1,2,\ldots,n-1$  

となる。

この差分方程式を求める方法を差分法(Finite-Difference Method)という。このとき,局所打切り誤差が $ O(h^{2} + k^{2})$であるように $ u(x_{i},y_{j})$の近似$ w_{ij}$をとると, $ i = 1,2,\ldots,n-1$ $ j = 1,2,\ldots,m-1$に対して

$\displaystyle 2\left[\left(\frac{h}{k}\right)^{2} + 1\right] w_{ij} - (w_{i+1,j...
...}) - \left(\frac{h}{k}\right)^{2}(w_{i,j+1} + w_{i,j-1}) = -h^{2}f(x_{i},y_{j})$ (53)

ただし,
$\displaystyle w_{0j}$ $\displaystyle =$ $\displaystyle g(x_{0},y_{j}), \ w_{nj} = g(x_{n},y_{j}), \ j = 0,1,\ldots,m$  
$\displaystyle w_{i0}$ $\displaystyle =$ $\displaystyle g(x_{i},y_{0}), \ w_{im} = g(x_{i},y_{m}), \ i = 1,2,\ldots, n-1$  

式(5.7)では$ u(x,y)$の近似を一般に

$\displaystyle (x_{i-1},y_{j}), \ (x_{i},y_{j}),\ (x_{i+1},y_{j}), \ (x_{i},y_{j-1}), \ (x_{i},y_{j+1})$

の点で行う必要がある。これらの点を図に表すと

図 5.2:

\includegraphics[width=6cm]{NUMERICAL/figure12-5.eps}

となる。これより,これらの点は点 $ (x_{i},y_{j})$に対して星型となっている。ここで,境界条件を用いると,$ w_{ij}$についての $ (n-1)(m-1)$ $ (n-1)(m-1)$列の連立方程式を解くことになる。

この連立方程式は,内部メッシュ点を番号を振りなおすことによりもっと効率よく解くことができる。その方法として推薦されているのは,

$\displaystyle P_{l} = (x_{i},y_{j}), \ w_{l} = w_{i,j}$

ただし, $ l = i + (m-1-j)(n-1)$である。この方法によると,ラベルが左から右へ,そして上から下へとふられる。

図 5.3:

\includegraphics[width=5cm]{NUMERICAL/figure12-6.eps}

例題 5.1   隣合う2つの境界は0℃に保たれていて,残りの2つの境界は0℃から100℃まで線形に増加する境界条件のもとで,大きさ0.5m×0.5mの薄い鉄板の定常温度分布を求めよ。

$ x$軸と$ y$軸に沿って,零の境界を配置すると,この境界値問題は

$\displaystyle \frac{\partial^{2}u}{\partial x^{2}}(x,y) + \frac{\partial^{2}u}{\partial y^{2}}(x,y) = 0$

ここで, $ (x,y) \in R=\{(x,y) : 0 < x < 0.5, \ 0 < y < 0.5\}$で,境界条件は

$\displaystyle u(0,y) = 0, \ u(x,0) = 0, u(x,0.5) = 200x, \ u(0.5,y) = 200y$

と表せる。

ここで,$ n = m = 4$とおくと,この問題は次の図のような格子点を持ち,その差分方程式は,$ i = 1,2,3$$ j = 1,2,3$に対して

$\displaystyle 4w_{i,j} - w_{i+1,j} - w_{i-1,j} - w_{i,j-1} - w_{i,j+1} = 0$

で与えられる。

図 5.4:

\includegraphics[width=6cm]{NUMERICAL/figure12-7.eps}

これを新しくふられたラベルによる格子点 $ w_{i} = u(P_{i})$を用いて表すと,$ P_{i}$における式は,

\begin{displaymath}\begin{array}{lrl}
P_{1} : & 4w_{1} - w_{2} - w_{4} = & w_{0,...
...9} : & 4w_{9} - w_{8} - w_{6} = & w_{3,0} + w_{4,1}
\end{array}\end{displaymath}

となる。ここで,右辺は境界条件より

$\displaystyle w_{1,0} = w_{2,0} = w_{3,0} = w_{0,1} = w_{0,2} = w_{0,3} = 0$

$\displaystyle w_{1,4} = w_{4,1} = 25, \ w_{2,4} = w_{4,2} = 50,\ w_{3,4} = w_{4,3} = 75$

となるので,この問題から作られる連立方程式は

$\displaystyle \left(\begin{array}{rrrrrrrrr}
4 & -1 & 0 & -1 & 0 & 0 & 0 & 0 & ...
...}{c}
25\\
50\\
150\\
0\\
0\\
50\\
0\\
0\\
25
\end{array}\right)$

で表せる。これより $ w_{1},w_{2},\ldots,w_{9}$はGauss-Seidel法を用いて求めることができる。その結果は次のようになる。

i 1 2 3 4 5 6 7 8 9
$ w_{i}$ 18.75 37.50 56.25 12.50 25.00 37.50 6.25 12.50 18.75

Poisson方程式の差分法による近似アルゴリズム


Poisson方程式

$\displaystyle \frac{\partial^{2}u}{\partial x^{2}}(x,y) + \frac{\partial^{2}u}{\partial y^{2}}(x,y) = f(x,y), \ a \leq x \leq b, c \leq y \leq d$

の近似解を次の境界条件

$\displaystyle x = aまたはx = bでc \leq y \leq dのとき,u(x,y) = g(x,y) $

$\displaystyle y = cまたはy = dでa \leq x \leq bのとき,u(x,y) = g(x,y)$

で求める。

===============================================

\begin{displaymath}\begin{array}{l}
{\rm 入力}\hspace{2.6ex} {端点}\ a,b,c,d; {タ...
... 1,2,\ldots,m-1におけるu(x_{i},y_{j})の近似w_{i,j}.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(1)\hspace{4ex} h \leftarrow (b-a)/n;\\
\hspace{7ex} k = (d-c)/m;\\
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(2)\hspace{4ex} {\rm for}\ i =1,2,\ldots,n-1...
...rm do}\\
\hspace{11ex} \ x \leftarrow a + ih;\\
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(3)\hspace{4ex} {\rm for}\ j = 1,2,\ldots,m-...
...o}\\
\hspace{11ex} \ y_{j} \leftarrow c + jk.\\
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(4)\hspace{4ex} {\rm for}\ i = 1,2,\ldots,n-...
...{\rm do};\\
\hspace{15ex} \ w_{i,j} \leftarrow 0.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(5)\hspace{4ex} \ \lambda \leftarrow h^{2}/k...
... 2(1 + \lambda);\\
\hspace{7ex} l \leftarrow 1.\
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(6)\hspace{4ex} {\rm while}(l \leq N)\ {\rm ...
...-1}\vert;\\
\hspace{11ex} w_{1,m-1} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(8)\hspace{8ex} {\rm for}\ i = 2,3,\ldots,n-...
...z\vert;\\
\hspace{19ex} \ w_{i,m-1} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(9)\hspace{8ex} z \leftarrow (-h^{2}f(x_{n-1...
...vert;\\
\hspace{19ex} \ w_{n-1,m-1} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(10)\hspace{7ex} {\rm for}\ j = m-2,m-3,\ldo...
... W_{1,j} - z\vert;\\
\hspace{19ex} \ w_{1,j} = z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(12)\hspace{11ex} {\rm for}\ i = 2,3,\ldots,...
...} - z\vert;\\
\hspace{23ex} w_{i,j} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(13)\hspace{11ex} z \leftarrow (-h^{2}f(x_{n...
...- z\vert;\\
\hspace{19ex} w_{n-1,j} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(14)\hspace{7ex} z \leftarrow (-h^{2}f(x_{1}...
...- z\vert;\\
\hspace{15ex} \ w_{1,1} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(15)\hspace{7ex} {\rm for}\ i = 2,3,\ldots,n...
...} - z\vert;\\
\hspace{19ex} w_{i,1} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(16)\hspace{7ex} z \leftarrow (-h^{2}f(x_{i}...
...- z\vert;\\
\hspace{15ex} w_{n-1,1} \leftarrow z.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(17)\hspace{7ex} {\rm if}(NORM \leq TOL\ {\r...
...3ex} {\rm OUTPUT}.\\
(19)\hspace{7ex} {\rm STOP}.
\end{array}\end{displaymath}

\begin{displaymath}\begin{array}{l}
(20)\hspace{7ex} \ l \leftarrow l + 1.\\
(...
...辛ー鷽瑤鯆兇┐泙靴ソ');\\
\hspace{7ex} {\rm STOP}.
\end{array}\end{displaymath}

===============================================

例題 5.2   次のPoisson方程式の境界値問題の近似解を差分法を用いて求めよ。

$\displaystyle \frac{\partial^{2} u}{\partial x^{2}}(x,y) + \frac{\partial^{2} u}{\partial y^{2}}(x,y) = xe^{y}, \ 0 < x < 2,\ 0 < y < 1$


$\displaystyle u(0,y)$ $\displaystyle =$ $\displaystyle 0, u(2,y) = 2e^{y}, \ 0 \leq y \leq 1$  
$\displaystyle u(x,0)$ $\displaystyle =$ $\displaystyle x, u(x,1) = e^{x}, \ 0 \leq x \leq 2$  

上のアルゴリズムを用いて厳密解 $ u(x,y) = xe^{y}$$ n = 6$$ m = 5$で近似する。

i j $ x_{i}$ $ y_{j}$ $ w_{i,j}$
1 1 0.33333333 0.20000000 0.40725839
1 2 0.33333333 0.40000000 0.49747103
1 3 0.33333333 0.60000000 0.60758125
1 4 0.33333333 0.80000000 0.74199595
2 1 0.66666667 0.20000000 0.81451486
2 2 0.66666667 0.40000000 0.99494013
2 3 0.66666667 0.60000000 1.21516201
2 4 0.66666667 0.80000000 1.48399266
3 1 1.00000000 0.20000000 1.22175775
3 2 1.00000000 0.40000000 1.49238811
3 3 1.00000000 0.60000000 1.82272256
3 4 1.00000000 0.80000000 2.22597759
4 1 1.33333333 0.20000000 1.62895765
4 2 1.33333333 0.40000000 1.98976644
4 3 1.33333333 0.60000000 2.43021217
4 4 1.33333333 0.80000000 2.96791746
5 1 1.66666667 0.20000000 2.03603945
5 2 1.66666667 0.40000000 2.48695275
5 3 1.66666667 0.60000000 3.03749865
5 4 1.66666667 0.80000000 3.70971893

演習問題 5.1  


=2.6zw =1 1. 差分法を用いて次の楕円型偏微分方程式の近似解を求めよ。

$\displaystyle \frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} = 4,\ 0 \leq x \leq 1, 0 < y <2$

$\displaystyle u(x,0) = x^{2},\ u(x,2) = (x - 2)^{2},\ 0 \leq x \leq 1$

$\displaystyle u(0,y) = y^{2},\ u(1,y) = (y - 1)^{2},\ 0 \leq y \leq 2$

ただし,$ h = k = 1$とする。また,厳密解 $ u(x,y) = (x - y)^{2}$と比較せよ。

=2.6zw =1 2. 差分法を用いて次の楕円型偏微分方程式の近似解を求めよ。

$\displaystyle \frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} = 4,\ 0 \leq x \leq 1, 0 < y <2$

$\displaystyle u(x,0) = x^{2},\ u(x,2) = (x - 2)^{2},\ 0 \leq x \leq 1$

$\displaystyle u(0,y) = y^{2},\ u(1,y) = (y - 1)^{2},\ 0 \leq y \leq 2$

ただし,$ h = k = 1$とする。また,厳密解 $ u(x,y) = (x - y)^{2}$と比較せよ。


next up previous contents index
Next: 放物線型偏微分方程式 Up: 偏微分方程式 Previous: 3つの偏微分方程式   目次   索引
administrator 平成16年10月26日