変数分離法によるポアソン方程式の解法

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

問題設定および基礎式

2次元のポアソン方程式
\begin{align}
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = c
\end{align}の解を有限な矩形領域\(-a < x < a\)、\(-b < y < b\)で求めることを考える（ただし、\(c\)は定数）。境界条件は
\begin{align}
u(-a,y)=u(x,-b)=u(a,y)=u(x,b)=u_0
\end{align}とする。
この問題の解は以下のように表される。
\begin{align}
u(x,y)=u_p(x,y)+u_c(x,y)+u_0
\end{align}ここで、\( u_p \)は与式および同次境界条件を満たす特殊解、\( u_c \)は与式の同次方程式を満たす解である。ここで、同次境界条件とは
\begin{align}
u(-a,y)=u(x,-b)=u(a,y)=u(x,b)=0
\end{align}であり、ポアソン方程式の同次方程式とはラプラス方程式
\begin{align}
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0
\end{align}である。

変数分離法による求解

まず、\( u_c \)を変数分離法によって求める。\( u_c(x,y)=X(x)Y(y) \)を解として仮定し、ラプラス方程式に代入すると、次式を得る。
\begin{align}
\frac{1}{X}\frac{d^2 X}{d x^2} = - \frac{1}{Y}\frac{d^2 Y}{d y^2} = k
\end{align}ここで、\( k \)は分離定数であり、今の場合、\( k \neq 0 \)である。そこで、\( k=\lambda^2 \)とおくと、次の2つの常微分方程式を得る。
\begin{align}
X'' - \lambda X &= 0 \\
Y'' + \lambda Y &= 0
\end{align}これらの解はただちに、
\begin{align}
X &= C_1 e^{\lambda x} + C_2 e^{-\lambda x} \\
Y &= D_1 \cos(\lambda y) + D_2 \sin(\lambda y)
\end{align}となる。\( Y \)に関する境界条件から\( D_1 \neq 0 \)とおくと、
\[ \lambda_n = \frac{(2n-1)\pi}{2b} \]となり、これより\( D_2=0 \)。\( X \)に関する境界条件から、
\[ C_1+C_2=0 \ \lor \ C_1-C_2=0 \]となるので、\( C_1=C_2=C/2 \)とおくと、\( u=XY \)より、
\begin{align}
C D_1 \mathrm{cos h}\left( \frac{(2n-1)\pi x}{2b} \right) \cos\left( \frac{(2n-1)\pi y}{2b} \right)
\end{align}が解となる。ラプラス方程式は線形方程式であるため、一般解はこれらの解を重ね合わせたものになる。したがって、\( a_n \)を新たな任意定数として、
\begin{align}
u_c(x,y) = \sum_{n=1}^{\infty} a_n \mathrm{cos h}\left( \frac{(2n-1)\pi x}{2b} \right) \cos\left( \frac{(2n-1)\pi y}{2b} \right)
\end{align}
次に$a_n$を同次境界条件におけるポアソン方程式を満たすように決定する。次の関数
\begin{align}
\frac{c}{2}y^2+c_1 y+c_2
\end{align}は明らかにポアソン方程式の特殊解の1つである。同次境界条件より任意定数を定めると、
\begin{align}
u_p(x,y) = \frac{c}{2}(y^2 - b^2)
\end{align}となる。これより、$x$に関する境界条件から、
\begin{align*}
\begin{split}
u(\pm a,y) &= \frac{c}{2}(y^2 - b^2) + \sum_{n=1}^{\infty} a_n \mathrm{cos h}\left( \frac{(2n-1)\pi a}{2b} \right) \cos\left( \frac{(2n-1)\pi y}{2b} \right) \\
&= 0
\end{split}
\end{align*}上式の両辺に
\[ \cos\left( \frac{(2m-1)\pi y}{2b} \right) \]
をかけ、$-b < y < b$まで積分すると、三角関数の直交性より、
\begin{align}
\int_{-b}^{-b} \cos\left( \frac{(2n-1)\pi y}{2b} \right) \cos\left( \frac{(2m-1)\pi y}{2b} \right) dy = b \delta_{nm}
\end{align}ここで、$\delta_{nm}$はクロネッカーのデルタである。これより、
\begin{align}
a_m \mathrm{cos h}\left( \frac{(2m-1)\pi a}{2b} \right) &= - \frac{1}{b} \int_{-b}^{-b} \frac{c}{2}(y^2 - b^2) \cos\left( \frac{(2m-1)\pi y}{2b} \right) dy \\
&= - \frac{1}{b} \frac{c}{2} \frac{32b^3}{\pi^3} \frac{(-1)^m}{(2m-1)^3}
\end{align}したがって、求める解は$u(x,y)=u_p(x,y)+u_c(x,y)+u_0$より、
\begin{align}
u(x,y) &= \frac{c}{2}(y^2 - b^2) + u_0 \notag \\
&- \frac{16b^2 c}{\pi^3} \sum_{n=1}^{\infty} \frac{(-1)^n}{(2n-1)^3} \frac{ \mathrm{cos h}\left( \frac{(2n-1)\pi x}{2b} \right) }{ \mathrm{cos h}\left( \frac{(2n-1)\pi a}{2b} \right) } \cos\left( \frac{(2n-1)\pi y}{2b} \right)
\end{align}

計算例

$c=-1$、$a=b=0.5$、$u_0=0$として、$n=50$までの和を計算したときの解の3Dプロットを以下に示す。

グラフの作成にはgnuplotを使用した。以下にコードを示す（級数部分の計算はここを参考にした）。

# 級数部分の定義
series(x,y,n) = sum[i=1:n] (-1)**(i)/(2*i-1)**3*cosh((2*i-1)*pi*x/(2*b))/cosh((2*i-1)*pi*a/(2*b))*cos((2*i-1)*pi*y/(2*b))
# ポアソン方程式の解
u(x,y,n)=c/2*(y**2-b**2)+u0-16*b**2*c/pi**3*series(x,y,n)

c=-1.0
a=0.5
b=0.5
u0=0.0

set xlabel 'x' font "Times-Italic"
set ylabel 'y' font "Times-Italic"
set xrange [-a:a]
set yrange [-b:b]
set pm3d at b
set view 60,30
set isosamples 100
set palette rgbformulae 22,13,-31

set terminal postscript enhanced color
set output "sol3D.eps"
splot u(x,y,50) notitle
set output

まとめ

上のグラフを見ると明らかに$y$軸上での境界条件が満たされていない（メッシュが曲がっている）。Mathematicaで計算したときには同じ項数でも精度よく境界条件が満たされていた（メッシュもまっすぐだった）ので、これはgnuplotに特有の問題かもしれない（原因？　知るか！）。最初はMathematicaの計算結果をTableで書き出してやろうかと思ったのだが、gnuplotで級数の計算ができると聞き、どうせグラフはgnuplotで書くんだしと思ってやってみたのだが…なんだか微妙なオチになってしまった。答えとしては多分合っている…と思う（保証はない）。
まぁ元々、「ポアソン方程式　変数分離法」で検索しても具体的な計算例が見つからず、図書館で見つけた本を参考にこうして備忘録をネットに上げてやろうと思ったというのが動機なので、これでよしとする。ネットの情報の大半は静電ポテンシャルに対するグリーン関数の解説か、ラプラス方程式の変数分離法までしかない。結局のところ、ポアソン方程式の特殊解と合わせて任意定数を決定するという手順が一番知りたかったことだった。以下が参考書である（すばりそのままの例題が載ってて助かった）。

偏微分方程式 (マックォーリ 初歩から学ぶ数学大全)

追記（20140514）

グラフが境界条件を満たしていない理由がわかったので補足する。gnuplotは内部の計算に整数型と実数型の変数を区別して使っている。そのため、実数の計算に整数型の値を使うと結果がおかしくなる。上記のコードの級数の定義部分を

# 級数部分の定義
series(x,y,n) = sum[i=1:n] 1.0*(-1)**(i)/(2*i-1)**3*cosh((2*i-1)*pi*x/(2*b))/cosh((2*i-1)*pi*a/(2*b))*cos((2*i-1)*pi*y/(2*b))
# ポアソン方程式の解
u(x,y,n)=1.0*c/2*(y**2-b**2)+u0-16*b**2*c/pi**3*series(x,y,n)

と「1.0」をかけて実数型に変換することで正しい計算結果が得られる。もちろん、実数の数値に「16.」のようにピリオドをつけてもよい。以下が正しいグラフである。きちんと境界条件が満たされていることがわかる。

わかってみれば何とも初歩的なミスであった。というか、前にも同じようなことをしでかしたような気がする。しっかりしろよ、俺。