格子状に並んでいるわけでもないランダムなサンプル点に対して補間を行う手法の1つとして Radial Basis Function (RBF) 補間が知られています。その RBF には様々な応用があり、私の研究にも応用しているのですが、今回はその RBF を用いて偏微分方程式を数値的に解く手法について紹介します。

RBF 補間

最初に RBF 補間について簡単に紹介しておきます。

RBF 補間では、点 $r_1, r_2, \ldots, r_n$ におけるデータ $f_1, f_2, \ldots, f_n$ を補間するのに

f(r)=\sum_{i=1}^n c_i \phi(\|r-r_i\|_2)

のような関数を考えます。この式の $\phi$ が RBF と呼ばれる関数です。
$c_i$ は係数で次の方程式を解くことで求めます。

\begin{pmatrix}
\phi(\|r_1-r_1\|_2) & \phi(\|r_1-r_2\|_2) & \cdots & \phi(\|r_1-r_n\|_2) \\
\phi(\|r_2-r_1\|_2) & \phi(\|r_2-r_2\|_2) & \cdots & \phi(\|r_2-r_n\|_2) \\
\vdots & \vdots & \ddots & \vdots \\
\phi(\|r_n-r_1\|_2) & \phi(\|r_n-r_2\|_2) & \cdots & \phi(\|r_n-r_n\|_2)
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2 \\ \vdots \\ c_n
\end{pmatrix}
=
\begin{pmatrix}
f_1 \\ f_2 \\ \vdots \\ f_n
\end{pmatrix}

この方程式の係数行列は例えば

$\phi(r) = e^{-cr^2}$ (Gaussians)
$\phi(r) = (c^2 + r^2)^{\beta/2}$ (Multiquadrics)

のような RBF において正定値になることが知られています。つまり、方程式がただ1つの解を持ちます。ただし、一般に係数行列の条件は悪くなりやすいことが知られており、それへの対処法も色々と研究されています。

偏微分方程式への応用

ここでは簡単のため Laplace 方程式

\triangle u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0

を考え、境界上での $u$ が与えらえている状態（Dirichlet 境界条件）を考えます。
離散化のため境界の内側で $r_1, r_2, \ldots, r_n$ という点を取り、境界上で $r_{n+1}, r_{n+1}, \ldots, r_N$ という点を取り、次の式を考えます。

u(r) = \sum_{i=1}^n c_i \triangle \phi(\|r-r_i\|_2) + \sum_{i=n+1}^N c_i \phi(\|r-r_i\|_2)

ここで、通常の RBF 補間と比べると境界の内側の点に対しては Laplacian $\triangle$ がついています。これには係数 $c_i$ を求める方程式の対称性をよくする意味があります。実際、$c_i$ を求める方程式は次のように書けます。

\begin{pmatrix}
\triangle\triangle\phi(\|r_1-r_1\|_2)  & \cdots & \triangle\triangle\phi(\|r_1-r_n\|_2) &
\triangle\phi(\|r_1-r_{n+1}\|_2) & \cdots & \triangle\phi(\|r_1-r_N\|_2) \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\triangle\triangle\phi(\|r_n-r_1\|_2)  & \cdots & \triangle\triangle\phi(\|r_n-r_n\|_2) &
\triangle\phi(\|r_n-r_{n+1}\|_2) & \cdots & \triangle\phi(\|r_n-r_N\|_2) \\
\triangle\phi(\|r_{n+1}-r_1\|_2)  & \cdots & \triangle\phi(\|r_{n+1}-r_n\|_2) &
\phi(\|r_{n+1}-r_{n+1}\|_2) & \cdots & \phi(\|r_{n+1}-r_N\|_2) \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\triangle\phi(\|r_N-r_1\|_2)  & \cdots & \triangle\phi(\|r_N-r_n\|_2) &
\phi(\|r_N-r_{n+1}\|_2) & \cdots & \phi(\|r_N-r_N\|_2)
\end{pmatrix}
\begin{pmatrix}
c_1 \\ \vdots \\ c_n \\ c_{n+1} \\ \vdots \\ c_N
\end{pmatrix}
=
\begin{pmatrix}
0 \\ \vdots \\ 0 \\ u_{n+1} \\ \vdots \\ u_N
\end{pmatrix}

ここで、$u_{n+1}, \ldots, u_N$ は境界上の点 $r_{n+1}, \ldots, r_N$ における $u$ の値です。
また、$\triangle$ は中の変数 $r_i$ に作用させた状態、つまり

\triangle\phi(\|r_i-r_j\|_2) = \triangle\phi(\|r-r_j\|_2)|_{r=r_i}

を表すものとします。
この方程式の係数行列は $\phi(r) = e^{-cr^2}$ において正定値になることが証明されていますし、$(c^2 + r^2)^{\beta/2}$ においても $\beta > 4$ であれば正定値になるそうです。
そこで、$\phi(r) = e^{-cr^2}$ を使って計算をしてみました。
このとき、$\phi$の微分は次のようになります。

\begin{aligned}
\triangle\phi(\|r-r'\|_2)&=4c(c\|r-r'\|_2^2-1)e^{-c\|r-r'\|_2^2} \\
\triangle\triangle\phi(\|r-r'\|_2)&=16c^2(c^2\|r-r'\|_2^4-4c\|r-r'\|_2^2+2)e^{-c\|r-r'\|_2^2}
\end{aligned}

これを用いてプログラム（後ろに載せておきます）を作り、関数 $u(x,y) = xy$ を真の解としてテストしてみました。その結果が次の図になります。

True は正しい $u$ で、Calculated は RBF を用いて求めた $u$ で、Error は誤差です。
単純な実装ですが、3桁くらいの精度は出ているようです。

あとがき

今回は RBF を用いた偏微分方程式の解法について書きました。偏微分方程式を解く手法として有名なものでは有限要素法や境界要素法といったものもありますが、それらに比べて簡単に試せるのがこの手法の利点だと思います。気になった方やちょっと偏微分方程式を解きたくなった方は試してみてください。

参考文献

以下の文献はいずれも今回の内容を含んでいます。
個人的には2つ目の文献が一番おすすめです。

C. Franke and R. Schaback, "Solving Partial Difference Equations by Collocation using Radial Basis Functions," Applied Mathematics and Computations, Vol. 93, pp. 73-82, 1998.
G. E. Fasshauer, Meshfree Approximation Methods with MATLAB. World Scientific, 2007.
B. Fornberg and N. Flyer, A Primer on Radial Basis Functions with Applications to the Geosciences. SIAM, 2015.

ソースコード

以下に本文中の実験で用いた Matlab のソースコードを載せておきます。
方程式を解くのにわざわざ固有値分解を用いているのは、解きたい方程式の条件が悪いのを安定化するためです。

{Test_RBF_Laplace2d_20181209.m}

% test of solving 2D Laplace equation using RBF
clear variables; close all;
rng('default');

% boundary points
res = 10;
x = (0:(res-1))' / res;
pos_bound = [...
    x, zeros(res, 1); ...
    x, ones(res, 1); ...
    zeros(res, 1), x; ...
    ones(res, 1), x];

% inner points
pos_inner = rand(res^2, 2);

% true solution
u = @(pos) pos(:,1) .* pos(:,2);

% RBF
c = res^2 / 10;
dist2 = @(pos1, pos2) (pos1(:,1) - pos2(:,1)').^2 + (pos1(:,2) - pos2(:,2)').^2;
rbf = @(dist2) exp(-c * dist2);
rbf_dif1 = @(dist2) 4 * c * (c * dist2 - 1) .* exp(-c * dist2);
rbf_dif2 = @(dist2) 16 * c^2 * (c^2 * dist2.^2 - 4 * c * dist2 + 2) .* exp(-c * dist2);

% compute the coefficients for RBF
K_ii = rbf_dif2(dist2(pos_inner, pos_inner));
K_ib = rbf_dif1(dist2(pos_inner, pos_bound));
K_bb = rbf(dist2(pos_bound, pos_bound));
K = [K_ii, K_ib; K_ib', K_bb];
f = [zeros(res^2, 1); u(pos_bound)];
[V, d] = eig(K, 'vector');
coeff = V * diag(d.^(-1)) * V' * f;

% calculate the results on test points and plot them
res = 50;
x = (0:(res-1) + 0.5) / res;
y = x;
[X, Y] = meshgrid(x, y);
pos_test = [reshape(X, res^2, 1), reshape(Y, res^2, 1)];
u_test = [rbf_dif1(dist2(pos_test, pos_inner)), rbf(dist2(pos_test, pos_bound))] ...
    * coeff;
u_true = u(pos_test);
U_test = reshape(u_test, res, res);
U_true = reshape(u_true, res, res);
fprintf('Relative Error: %e\n', ...
    norm(U_test - U_true, 'fro') / norm(U_true, 'fro'));

h = figure;
temp = h.Position;
temp(1) = temp(1) * 0.3;
temp(3) = temp(3) * 2;
h.Position = temp;

subplot(1, 3, 1);
imagesc(x, y, U_true);
caxis([0 1]);
colorbar;
axis equal tight xy;
xlabel('x'); ylabel('y');
title('True');

subplot(1, 3, 2);
imagesc(x, y, U_test);
caxis([0 1]);
colorbar;
axis equal tight xy;
xlabel('x'); ylabel('y');
title('Calculated');

subplot(1, 3, 3);
imagesc(x, y, abs(U_true - U_test));
colorbar;
axis equal tight xy;
xlabel('x'); ylabel('y');
title('Error');