%=====================================================================
% 秋田高専 5E 計算機応用テキスト
%　　テーマ　偏微分方程式(ラプラス方程式)
%    本テキストの内容
%
% last updated 2005.2.7
%    created by  Masashi Yamamoto
%     e-mail yamamoto@akita-nct.jp
%=====================================================================
\documentclass[10pt,a4paper]{jarticle}
\usepackage{html,graphicx,amsmath,amssymb,ascmac,float}
\newcommand{\vm}[1]{\boldsymbol{#1}}
\newcommand{\hugesymbol}[1]{\mbox{\strut\rlap{\smash{\Huge$#1$}}\quad}}
\newcommand{\tw}[1]{\texttt{#1}}
\oddsidemargin 0mm  %左の余白 25.4mm-0mm　奇数ページ
\evensidemargin 0mm %左の余白 25.4mm-0mm　偶数ページ
\textwidth 160mm
%
\begin{document}
\title{偏微分方程式(ラプラス方程式)}
\author{山本昌志\thanks{独立行政法人　秋田工業高等専門学校　電気工学科}}
\date{2006年2月3日}
\maketitle
%
%
%=====================================================================
\section{はじめに}
%=====================================================================
以前，常微分方程式の数値計算について学習した．独立変数が1個のものを常微分方程式，
2個以上のものを偏微分方程式と言うのは数学の授業で学んだとおりである．実際，自然
現象は常微分方程式よりも，偏微分方程式で記述されることが多い．常微分方程式が役に
立たないと言っているのではなく，より広範囲には偏微分方程式が使われているというこ
とである．自然界が，$(x, y, z)$の3次元と時間$t$を合わせた4次元で成り立っているた
めである．

ここでは，偏微分方程式，特にラプラス方程式を差分法というテクニックで数
値計算する方法を学習する．偏微分方程式は，いろいろなものがあるが，最初
に学習する分には，意味がわかりやすい方程式と言うことで，これを教
材に選んだ．実際には，図\ref{fig:potential}の静電磁場や，図
\ref{fig:temperature}の熱の問題に，この方程式は表れる．
%
\begin{itemize}
 \item 図\ref{fig:potential}は，紙面と垂直方向に無限に長い正方形の金属筒に，電線
       が2本通っている．正方形の筒は0Vにアースされており，2本の電線はそれぞれ，
       30Vと-20Vである．この状態で，筒内部のポテンシャル(電圧)の分布を求めなさい
       という問題である．
 \item 図\ref{fig:temperature}は，紙面方向に非常に長い豆腐があり，その周りは0
       ℃の水で満たされている．そして，その豆腐に30℃と-20℃の金属棒が突き刺
       さっている状況である．豆腐内部の温度分布を求めなさいというのが問題であ
       る．
\end{itemize}


これらは，物理的には異なる問題であるが，ポテンシャルや温度が満たす方程式は同じで
ある．方程式が同じならば解は同じで，同じ計算手法が使えることを理解して欲しい．こ
れらが満たすのはラプラス方程式
%
\begin{align}
 \nabla^2 \phi = 0
\end{align}
%
と言われるものである．$(x,y,z)$の3変数の偏微分方程式である．ここでは，
3次元の問題は大変なので，図\ref{fig:potential}や\ref{fig:temperature}のように紙
面の方向($z$方向)には一様とする．そのため，$z$方向の微分はゼロとなるので，ラプラ
ス方程式は，
%
\begin{align}
 \frac{\partial^2 \phi}{\partial x^2}+
 \frac{\partial^2 \phi}{\partial y^2}=0
 \label{eq:2D_Laplace}
\end{align}
%
と2次元問題になる．ここではこの偏微分方程式の近似解を数値計算により求めるのが目
的である．ここでの学習を通して，プログラムが完成すると，図\ref{fig:solution}のよ
うな解のグラフを求めることができる．
%
\begin{figure}[hbtp]
 \begin{tabular}{cc}
  \begin{minipage}[b]{0.45\hsize}
   \begin{center}
    \includegraphics[keepaspectratio, scale=0.8]{figure/potential.eps}
    \caption{ポテンシャルを求める問題．}
    \label{fig:potential}
   \end{center}
  \end{minipage} &
  \begin{minipage}[b]{0.45\hsize}
   \begin{center}
    \includegraphics[keepaspectratio, scale=0.8]{figure/temperature.eps}
    \caption{温度を求める問題．}
    \label{fig:temperature}
   \end{center}
  \end{minipage}
 \end{tabular}
\end{figure}
%
%
\begin{figure}[hbtp]
 \begin{center}
  \includegraphics[keepaspectratio, scale=1.0]{figure/solution.eps}
  \caption{差分法により計算された，ポテンシャルや温度のグラフ．}
  \label{fig:solution}
 \end{center}
\end{figure}
%
%=====================================================================
\section{差分法による偏微分方程式の数値計算}
%=====================================================================
%---------------------------------------------------------------------
\subsection{差分方程式}
%---------------------------------------------------------------------
2次元のラプラス方程式を数値計算で近似解を求めることを考える．まずは，いつものよう
に，解$\phi(x,y)$をテイラー展開する．xおよび，y方向に微小変位$\pm h$があっ
た場合，
%
\begin{align}
 \phi(x+h,y)&=\phi(x,y)
 +\frac{\partial\phi}{\partial x}h
 +\frac{1}{2!}\frac{\partial^2\phi}{\partial x^2}h^2
 +\frac{1}{3!}\frac{\partial^3\phi}{\partial x^3}h^3
 +\frac{1}{4!}\frac{\partial^4\phi}{\partial x^4}h^4
 +\cdots\\
 \phi(x-h,y)&=\phi(x,y)
 -\frac{\partial\phi}{\partial x}h
 +\frac{1}{2!}\frac{\partial^2\phi}{\partial x^2}h^2
 -\frac{1}{3!}\frac{\partial^3\phi}{\partial x^3}h^3
 +\frac{1}{4!}\frac{\partial^4\phi}{\partial x^4}h^4
 -\cdots
\end{align}
%
となる．これらの式の辺々を足し合わせえると，
%
\begin{align}
 \left.\frac{\partial^2\phi}{\partial x^2}\right|_{x,y}
 &=\frac{1}{h^2}\left[
 \phi(x+h,y)-2\phi(x,y)+\phi(x-h,y)\right]-O(h^2)
 \label{eq:approx_d2fdx2}
\end{align}
%
が得られる．このことから，2階の偏導関数の値は微小変位$h$の場所の関数の値を用いて，
$h^2$の精度で近似計算ができることが分かる．すなわち，式(\ref{eq:approx_d2fdx2})
の右辺の第1項を計算すればよいのである．同じことをy方向についても行うと
%
\begin{align}
 \left.\frac{\partial^2\phi}{\partial y^2}\right|_{x,y}
 &=\frac{1}{h^2}\left[
 \phi(x,y+h)-2\phi(x,y)+\phi(x,y-h)\right]-O(h^2)
 \label{eq:approx_d2fdy2}
\end{align}
%
が得られる．

これらの式(\ref{eq:approx_d2fdx2})と(\ref{eq:approx_d2fdy2})を元の2次元ラプラス
方程式(\ref{eq:2D_Laplace})に代入すれば，
%
\begin{align}
 \phi(x+h,y)+\phi(x-h,y)+\phi(x,y+h)+\phi(x,y-h)-4\phi(x,y)=0
 \label{eq:sabun}
\end{align}
%
となる．これが，2次元ラプラス方程式の差分の式である．この式を眺めると，座標
$(x,y)$のポテンシャルの値$\phi(x,y)$は，周りの値の平均であることがわかる．

実際にこの式を数値計算する場合，例えば図\ref{fig:potential}のポテンシャルを求め
る時には，図\ref{fig:reagine_lattice}のように格子状\footnote {この格子のことをメッ
シュ(mesh)と言う事もある．}に区切り，その交点での値を求めることになる．ここでは，
xおよびy方向には等間隔$h$ で区切り計算を進めるが，等間隔である必要はない．多少，
式(\ref{eq:sabun}) は異なるが同じような計算は可能である．これまでの説明が理解で
きていれば，xとy方向の間隔が異なっても，式(\ref{eq:sabun})に対応する差分の式が作
れるはずである．
%
\begin{figure}[hbtp]
 \begin{center}
  \includegraphics[keepaspectratio, scale=1.0]{figure/lattice.eps}
  \caption{解くべき領域を格子に分割}
  \label{fig:reagine_lattice_all}
 \end{center}
\end{figure}
%

数値計算をする場合，$\phi(x,y)$や$\phi(x+h,y)$の形は不便なので，形式を
改める．図\ref{fig:reagine_lattice_all}の左下の座標を(0,0)として，格子点で
のポテンシャルを
%
\begin{equation}
 \begin{aligned}
  \phi(x,y)&=\phi(ih,jh)\\
  &=U_{i\,j}
 \end{aligned}
\end{equation}
%
とする．このようにすると，式(\ref{eq:sabun})は
%
\begin{align}
 U_{i+1\,j}+U_{i-1\,j}+U_{i\,j+1}+U_{i\,j-1}-4U_{i\,j}=0
 \label{eq:sabun_U}
\end{align}
%
となり，数値計算し易い形になる．このようにした場合の各格子点の様子を図
\ref{fig:reagine_lattice}に示す．

次の節で述べる境界条件を考えないとすると，ラプラス方程式は式(\ref{eq:sabun_U})の
連立方程式を解くだけである．格子に領域を分割することにより，難しげな偏微分方程式
が連立方程式に還元されたわけである．

%
\begin{figure}[hbtp]
 \begin{center}
  \includegraphics[keepaspectratio, scale=1.0]{figure/sabun.eps}
  \caption{差分の格子}
  \label{fig:reagine_lattice}
 \end{center}
\end{figure}
%
%---------------------------------------------------------------------
\subsection{境界条件(外周)}
%---------------------------------------------------------------------
実際に，連立方程式(\ref{eq:sabun_U})を計算する場合，困った問題が生じる．このまま
だと，式の数と未知数の数が合わないのである．たとえば，図\ref{fig:eq_at_boundary}
に示す境界を考える．すると，境界が式\ref{eq:sabun_U}の$(i,\,j)$になる場合，式が
作れないのである．すると，図\ref{fig:reagine_lattice}の中の電極が無い場合，可能
な連立方程式の数は，$(N_x-1)\times(N_y-1)$個である．未知数の数は，
$(N_x+1)\times(N_y+1)$個である．未知数の数の方が，$2(N_x+N_y)$個多いのである．そ
のため，予め，この余分の未知数となっている値を決めなくてはならない．実際これは，
偏微分方程式の境界条件を決めることに相当する．
%
\begin{figure}[hbtp]
 \begin{center}
  \includegraphics[keepaspectratio, scale=1.0]{figure/eq_boundary.eps}
  \caption{境界付近の差分方程式の可能性}
  \label{fig:eq_at_boundary}
 \end{center}
\end{figure}
%

そこで，境界上の格子点$2(N_x+N_y)$個の値を予め決める．こうすれば，式の数を減らさ
ないで，未知数の数を減らすことができる．要するに偏微分方程式を解くときの境界条件
を決めるのと同じ．驚いたことに，外周部(境界条件)の格子点の数が，ちょうど，不足し
ている方程式の数と同じなのである．自然は，都合良くできているのである．

懸命な諸君であれば，予め決める値は外周の境界上の格子点でなくても良いと考えるだろ
う．しかし，内部の点の値を決めてしまうと，連立方程式が1個減ってしまうので，未知
数と式の数の差は変わらない．これについては，良い説明が思い浮かばなかったので，そ
ういうものだと思ってください．
%
%---------------------------------------------------------------------
\subsection{境界条件(内部の電極)}
%---------------------------------------------------------------------
先ほどの説明通り内部の格子点のポテンシャルを決めてしまうと，その数だけ方程式が減
少する．したがって，必要なだけ内部のポテンシャルを決めても，式と未知数の数は同
じで，連立方程式は解ける．そのような理由で，図\ref{fig:reagine_lattice}
の電極内部の格子点のポテンシャル$U_{i\,j}$ は，問題で与えられたとおり，30Vと-20V
と固定しても，連立方程式の数は変わらない．
%
%=====================================================================
\section{実際の計算方法}
%=====================================================================
%---------------------------------------------------------------------
\subsection{ガウス・ザイデル法}
%---------------------------------------------------------------------
境界条件，境界でのポテンシャルの値を決めることで，方程式と未知数の数が一致するこ
とは前に示したとおりである．したがって，連立方程式を解けば，格子点のポテンシャル
は分かることになる．しかし，実際問題この方程式を作るのが大変である．未知数のポテ
ンシャルは，分かりやすくするために$U_{i\,j}$と行列で表示したが，実際に連立方程式
を解く場合，未知数としてベクトルになる．この式を書くのは大変厄介である．本当に大
変かどうかは，諸君が実際に，係数行列と同次項を求めてみれば分かる．

そこで，次のお手軽な方法で連立方程式を解くことにする．この方法は，連立方程式の係
数行列や同次項を求める必要がなく，プログラムは非常に簡単になる．
\begin{enumerate}
\item{外周の境界上の格子点のポテンシャルの値を$U_{i\,j}$に代入する．こ
     の値は，これ以降，ずっと一定とする．}
\item{電極の内部の格子点のポテンシャルの値を$U_{i\,j}$を代入する．この
     値も，これ以降，ずっと一定とする．}
\item{計算すべき内部の格子点のポテンシャルを，
     %
     \begin{align}
      U_{i\,j}=\frac{1}{4}
      \left[U_{i+1\,j}+U_{i-1\,j}+U_{i\,j+1}+U_{i\,j-1}\right]
      \label{eq:sabun_numerical_cal_U}
     \end{align}
     %    
     に従い計算する．}
\item{$U_{i\,j}$が収束するまで，繰り返し，式
     (\ref{eq:sabun_numerical_cal_U})を計算する．}
\end{enumerate}

実際，この方法で無限の回数ループ処理をすれば，真の解に収束するはずである．このよ
うに，反復により解に収束させる方法を\htmladdnormallink{反復法}
{http://www.akita-nct.jp/~yamamoto/lecture/2005/5E/Linear_eauations/relaxation_html/relaxation.html}
と言う．とくに，ここで示した方法を\htmladdnormallink{ガウス・ザイデル法}
{http://www.akita-nct.jp/~yamamoto/lecture/2005/5E/Linear_eauations/relaxation_html/node5.html}
という．以前学習したしたように，\htmladdnormallink{SOR法}
{http://www.akita-nct.jp/~yamamoto/lecture/2005/5E/Linear_eauations/relaxation_html/node6.html}
の方が収束が早いが，ここではプログラムが簡単なガウス・ザイデル法で計算するのが良
いであろう．
%
%---------------------------------------------------------------------
\subsection{内部電極の決め方}
%---------------------------------------------------------------------
ある注目している格子点が，ポテンシャルが固定されている内部電極の内側にあるか否か
判断しなくてはならない．内部にあるとそれは固定点で，先に示したガウス・ザイデル法
で値を変えてはならない．一方，外部だと境界でない限り，ガウス・ザイデル方でしゅう
そくするまで，計算を繰り返すことになる．そのようなことから，電極内部にあるか否か
は予め判断する必要がある．その判断は簡単である．格子点と電極中心の距離と電極半径
を比べれば良いのである．
%
%---------------------------------------------------------------------
\subsection{固定及び境界点の計算を省く方法}
%---------------------------------------------------------------------
固定電極内部の点や外部の境界の格子点では，ガウスザイデル方で計算する必要がない．
それを計算しないようにプログラムを書かなくてはならない．境界の計算を省くのは簡単
である．境界の格子点のポテンシャルは，$i=0, i=N_x, j=0, j=N_y$の場合である．この
場合，ガウス・ザイデル法で式(\ref{eq:sabun_numerical_cal_U})の$U_{i\,j}$を計算し
なければ良いのである．

次に，電極内部の点であるが，これは予め目印を付ければ良い．例えば，整数型の配列を
用意して，その値が1の場合は，電極の内部と目印をしておく．1以外の時，式
(\ref{eq:sabun_numerical_cal_U})の$U_{i\,j}$を計算する．

以上をまとめると，連立方程式を解く場合には，次のようなループとすれば良いであろう．
%
\begin{quote}
 \setlength{\baselineskip}{12pt}
 \begin{verbatim}
	for(k=1; k<=ガウス・ザイデル法の計算回数;k++){
	  for(j=1; j<= x方向の分割数-1 ; j++){
	    for(i=1; i<= y方向の分割数-1 ;i++){
	      if(フラグ[i][j] != 1){
	        u[i][j]=0.25*(u[i][j-1]+u[i][j+1]+u[i-1][j]+u[i+1][j]);
	      }
	    }
	  }
	}
 \end{verbatim}
\end{quote}
%

このループを図\ref{fig:cal_points}に示す．先のループがどのように計算されるか，よ
く考えよ．
%
\begin{figure}[hbtp]
 \begin{center}
  \includegraphics[keepaspectratio, scale=1.0]{figure/cal_points.eps}
  \caption{計算する格子点と順序}
  \label{fig:cal_points}
 \end{center}
\end{figure}
%
%
%=====================================================================
\section{練習問題}
%=====================================================================
%---------------------------------------------------------------------
\subsection{ポテンシャルの計算}
%---------------------------------------------------------------------
図\ref{fig:potential}のポテンシャル分布を求めるプログラムを作成せよ．ただし，詳
細な寸法は，図\ref{fig:question}の通りとする．
%
\begin{figure}[hbtp]
 \begin{center}
  \includegraphics[keepaspectratio, scale=1.0]{figure/question.eps}
  \caption{練習問題の境界条件とポテンシャル}
  \label{fig:question}
 \end{center}
\end{figure}
%
%---------------------------------------------------------------------
\subsection{ヒント}
%---------------------------------------------------------------------
途中まで作成したソースプログラムを私の講義ノートのWEBページに置いておく．プログラム
作成の手がかりがつかめない者は，これを完成させよ．

このヒントのプログラムの変数は，次のように使われている．
%
\begin{quote}
 \begin{description}
  \item[\tw{u[i][j]}]\quad\tw{i,j}番目の格子点のポテンシャルを入れる配列
  \item[\tw{x[i][j],y[i][j]}]\quad\tw{i,j}番目の格子点の座標を入れる配列
  \item[\tw{flag[i][j]}]\quad\tw{i,j}番目の格子点が，内部電極内ならば1，電極外な
	     らば1となるフラグ(目印)
  \item[\tw{iteration}]\quad ガウス・ザイデル法での最大計算回数
  \item[\tw{nlat}]\quad x方向とy方向の分割数．問題が正方形領域なので，同一としている．
 \end{description}
\end{quote}
%

また，使われている関数の動作は以下の通りである．
%
\begin{quote}
 \begin{description}
  \item[\tw{initialize}]\quad ポテンシャルとフラグの値をゼロに初期化する．
  \item[\tw{set\_cordinate}]\quad 格子点の座標を設定する．
  \item[\tw{set\_boundary\_wall\_pot}]\quad 外部境界の格子点のポテンシャルを
	     0，フラグを1に設定する．
  \item[\tw{set\_circle}]\quad 電極内部の格子点のポテンシャルを設定する．さらにその
	     フラグを1に設定する．
  \item[\tw{plot\_3d}]\quad 結果を3次元グラフに表示する．
 \end{description}
\end{quote}
%=====================================================================
%
\end{document}