代用電荷法（基本解近似解法）

代用電荷法（基本解近似解法）
Charge Simulation Method (Method of Fundamental Solutions)

periodic Stokes flow

２次元正方形格子上に並んだ球列をすぎるストークス流（基本解近似解法により計算）。

Stokes flow past a planar square array of spheres computed by the fundamental solution method based on the periodic fundamental solution.I am grateful to Mr. Daisuke Matsuoka (Graduate School of Science and Engineering, Ehime University) for his help in making this figure.

科学技術問題を数理的に扱おうとすると、必ずといっていいほど、偏微分方程式があらわれます。偏微分方程式の解を厳密に得る場合はまれで、ほとんどの場合、数値的に解を求めることになります。偏微分方程式の数値解法としては、有限要素法、差分法が代表的なもので多くの研究がなされています。私達の研究室で研究している代用電荷法も偏微分方程式の数値解法のひとつであり、とくにポテンシャル問題（流体力学、電磁気学などに現れます）によく用いられます。そして、次のようなよい性質を持っています。

プログラミングが簡単。
計算量が少ない。
ある条件下では真の解への収束が速い。

ポテンシャル問題に対する代用電荷法

では、これから代用電荷法について具体的に説明します。 $\mathcal{D}$ を2次元平面または3次元空間内の領域として、領域 $\mathcal{D}$ における次のポテンシャル問題、すなわち、 ラプラス方程式の境界値問題を考えます。

$\begin{cases} \Delta u = 0 & \quad \mbox{in} \ \mathcal{D}\\ u = f & \quad \mbox{on} \ \partial\mathcal{D} \end{cases}$

ここで、 $\partial\mathcal{D}$ は領域 $\mathcal{D}$ の境界、 $f$ は境界 $\partial\mathcal{D}$ 上で与えられた関数です。代用電荷法では、この境界値問題の解を次の形で近似します。

$u(\boldsymbol{x})\simeq u_N(\boldsymbol{x}) = \begin{cases} \displaystyle -\frac{1}{2\pi}\sum_{j=1}^{N}Q_j\log\|\boldsymbol{x}-\boldsymbol{\xi}_j\| & \mbox{(two-dimensional case)} \\ \displaystyle \frac{1}{4\pi}\sum_{j=1}^{N}\frac{Q_j}{\|\boldsymbol{x}-\boldsymbol{\xi}_j\|} & \mbox{(three-dimensional case)} \end{cases}$

csm

図1：代用電荷法のイメージ

ここで、 $\boldsymbol{\xi}_1,\ldots,\boldsymbol{\xi}_N$ は領域 $\mathcal{D}$ にユーザによって与えられた点で、 電荷点と呼びます。 $Q_1,\ldots,Q_N$ は実数係数で、電荷と呼びます。電荷の定め方についてはすぐ後で説明します。物理の言葉で言えば、近似解 $u_N(\boldsymbol{x})$ は点 $\boldsymbol{\xi}_j \ ( \: j = 1, \ldots, N \: )$ に置かれた点電荷 $Q_j$ がつくるクーロン電位の重ね合わせによって得られる電位に相当します。このように、代用電荷法は物理的直感に基づいた数値解法です（図1参照）。

ここで、近似解 $u_N(\boldsymbol{x})$ は領域 $\mathcal{D}$ でラプラス方程式 $\triangle u_N = 0$ を満たすことに注意しましょう。つまり、代用電荷法においては、近似解は問題となっている偏微分方程式（ここではラプラス方程式）を厳密に満たすのです。一方、境界条件については、条件を少し緩めて、次の拘束条件を満たすようにします。

$u_N(\boldsymbol{x}_i) = f(\boldsymbol{x}_i) \quad ( \: i = 1, \ldots, N \: ),$

ここで、 $\boldsymbol{x}_1,\ldots,\boldsymbol{x}_N$ は境界 $\partial\mathcal{D}$ 上にユーザによって与えられた点であり、 拘束点と呼びます（表1参照）。拘束条件は電荷 $Q_j$ に対する連立1次方程式

$\begin{bmatrix} G_{11} & \cdots & G_{1N} \\ \vdots & & \vdots \\ G_{N1} & \cdots & G_{NN} \end{bmatrix} \begin{bmatrix} Q_1 \\ \vdots \\ Q_N \end{bmatrix} = \begin{bmatrix} f(\boldsymbol{x}_1) \\ \vdots \\ f(\boldsymbol{x}_N) \end{bmatrix},$ $G_{ij} = \begin{cases} \displaystyle - \frac{1}{2\pi}\log\|\boldsymbol{x}_i - \boldsymbol{\xi}_j\| & \mbox{(two-dimensional case)} \\ \displaystyle \frac{1}{4\pi}\frac{1}{\|\boldsymbol{x}_i - \boldsymbol{\xi}_j\|} & \mbox{(three-dimensional case)} \end{cases} \quad ( \: 1 \leq i, j \leq N \: )$

を成し、この連立1次方程式を解くことにより電荷 $Q_j$ が定まります。

代用電荷法の一般化

これまで説明したポテンシャル問題に対する代用電荷法は，他の問題（波動問題，流体力学，弾性問題など）に応用することができます。

ポテンシャル問題に対する代用電荷法では，近似解は

$u_N(\boldsymbol{x}) = \sum_{j=1}^{N} Q_j E(\boldsymbol{x}-\boldsymbol{\xi}_j), \quad E(\boldsymbol{x}) = \begin{cases} \displaystyle - \frac{1}{2\pi}\log\|\boldsymbol{x}\| & \mbox{(two-dimensional case)} \\ \displaystyle \frac{1}{4\pi}\frac{1}{\|\boldsymbol{x}\|} & \mbox{(three-dimensional case)} \end{cases}$

と表され， $E(\boldsymbol{x})$ が $\boldsymbol{x}\neq\boldsymbol{0}$ で $\triangle E(\boldsymbol{x}) = 0$ を満たす，もっと正確に言えば，

$\triangle E(\boldsymbol{x}) = \delta(\boldsymbol{x})$ 　　　　（ $\delta(\boldsymbol{x})$ はデルタ関数）

を満たすことから，近似解 $u_N(\boldsymbol{x})$ がラプラス方程式 $\triangle u_N = 0$ を厳密に満たします。

この点に着目して，代用電荷法の一般化を考えてみます。 $\mathcal{L}$ を線形偏微分作用素， $\mathcal{D}$ を2次元平面または3次元空間内の領域として，次の斉次偏微分方程式の境界値問題を考えます。

$\begin{cases} \mathcal{L}u = 0 & \mbox{in} \ \mathcal{D} \\ u = f & \mbox{on} \ \partial\mathcal{D}. \end{cases}$

この偏微分方程式に対し， $E(\boldsymbol{x})$ を方程式

$\mathcal{L}E(\boldsymbol{x})=\delta(\boldsymbol{x})$

の解とします。 $E(\boldsymbol{x})$ を偏微分作用素 $\mathcal{L}$ の基本解と呼びます。そして、いま考えている偏微分方程式の境界値問題に対し、次の形の近似解を考えます。

$u(\boldsymbol{x}) \simeq u_N(\boldsymbol{x}) = \sum_{j=1}^{N} Q_j E(\boldsymbol{x} - \boldsymbol{\xi}_j).$

ここで， $\boldsymbol{\xi}_1, \ldots, \boldsymbol{\xi}_N$ は領域 $\mathcal{D}$ は外部にユーザによって与えられた点であり， ソース点と呼びます。 $Q_1, \ldots, Q_N$ はすぐあとで述べる方法により定める定数係数です。

近似解 $u_N(\boldsymbol{x})$ は領域 $\mathcal{D}$ で偏微分方程式 $\mathcal{L}u_N = 0$ を厳密に満たすことに注意しましょう。
境界条件については、もとの条件を緩めた拘束条件

$u_N(\boldsymbol{x}_i) = f(\boldsymbol{x}_i) \quad ( \: i = 1, \ldots, N \: )$

を近似解に課すことにします。ここで、 $\boldsymbol{x}_1, \ldots, \boldsymbol{x}_N$ は境界 $\partial\mathcal{D}$ 上にユーザによって与えられた点であり，拘束点と呼びます。拘束条件は $Q_1,\ldots,Q_N$ に関する連立1次方程式

$G_{ij} = E(\boldsymbol{x}_i - \boldsymbol{\xi}_j) \quad ( \: 1\leq i, j\leq N \: )$ をなし、この連立1次方程式を解くことにより係数 $Q_1,\ldots,Q_N$ が定まります。

以上で述べたのが代用電荷法の一般化です。このように、ポテンシャル問題以外の問題に一般化した代用電荷法は基本解近似解法 (method of fundamental solutions)と呼ばれることが多いです。

例として、2次元波動問題に対する代用電荷法（基本解近似解法）の応用について述べます。
速度 $c$ で伝搬する2次元波は、2次元波動方程式

$\triangle U = \frac{1}{c^2}\frac{\partial^2 U}{\partial t^2},$ 　　　　ここで、　　　 $U = U(x, y, t), \quad \triangle = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ 　　　（ラプラシアン）

で記述されます。 $U(x, y, t) = \mathrm{e}^{\mathrm{i}\omega t}u(x, y)$ （ $\omega$ は角振動数）と変数分離すると、 ヘルムホルツ方程式

$(\triangle+k^2)u=0 \quad \left( \: k = \sqrt{\frac{\omega}{c}} \: \right)$

を得ます。波動問題はこのヘルムホルツ方程式を解く問題に帰着されることが多いです。そこで、 $\mathcal{D}$ を2次元平面内の外部領域として、次のヘルムホルツ方程式の境界値問題を考えます。

$\begin{cases} (\triangle + k^2)u = 0 & \mbox{in} \ \mathcal{D} \\ u = f & \mbox{on} \ \partial\mathcal{D} \\ \displaystyle \lim_{r\rightarrow\infty}\sqrt{r} \left(\frac{\partial u}{\partial r} - \mathrm{i}ku \right) = 0 & ( \: r = \|\boldsymbol{x}\| \: ), \end{cases}$

ここで、 $f$ は境界 $\partial\mathcal{D}$ 上で定義された関数、 3番目の式はゾンマーフェルトの外部放射条件であり、外部波動問題を考える時標準的に課される条件です。
いまの場合、ヘルムホルツ作用素 $\mathcal{L} = \triangle + k^2$ の基本解は、

$E(\boldsymbol{x}) = H_0^{(1)}(k\|\boldsymbol{x}\|)$ となります。 $H_0^{(1)}(\cdot)=J_0(\cdot)+\mathrm{i}Y_0(\cdot)$ は0次の第1種ハンケル関数です（ $J_0(\cdot)$ は0次の第1種ベッセル関数、 $Y_0(\cdot)$ は0次の第2種ベッセル関数（ノイマン関数）、図2参照）。

bessel

図2：ベッセル関数

したがって、上の境界値問題に対して代用電荷法（基本解近似解法）の近似解は次の形となります。

$u(\boldsymbol{x}) \simeq u_N(\boldsymbol{x}) = \sum_{j=1}^{N} Q_j H_0^{(1)}(k\|\boldsymbol{x} - \boldsymbol{\xi}_j\|),$

ここで、 $\boldsymbol{\xi}_1,\ldots,\boldsymbol{\xi}_N$ は領域 $\mathcal{D}$ 外部にユーザによって与えられたソース点、 $Q_1,\ldots,Q_N$ は定数係数であり、上に述べたように拘束条件を課すことにより定めます。
物理的なイメージでいうと、近似解 $u_N(\boldsymbol{x})$ は点 $\boldsymbol{\xi}_1,\ldots,\boldsymbol{\xi}_N$ に波源をもつ２次元球面波を重ねあわせたものとなります（図3参照）。

csm_wave

図3：波動問題に対する代用電荷法のイメージ

図4は、円柱によって平面波が散乱される様子を、基本解近似解法によって計算したものです（図では円柱は描かれていませんが）。

scattered_wave

図4：円柱による平面波の散乱

代用電荷法に関する我々の成果と現在の研究

緒方研究室における代用電荷法に関する研究業績には下記のものがあります。

空間周期的問題に対する代用電荷法の拡張。
上の図は、球体の３次元無限格子を過ぎる粘性流（ストークス流）の様子を表します。
圧縮性流問題に対する代用電荷法。
代用電荷法による流体解析（ポテンシャル流の圧力解析など）。
波動問題（ヘルムホルツ方程式）に対する代用電荷法，理論誤差評価と数値実験。
高速多重極展開およびGPUコンピューティングによる代用電荷法の高速計算。
代用電荷法による逆問題解析。