Network上の拡散方程式

はじめに

明治大学総合数理学部 Advent Calender 2017の6日目の記事。やっとこさ人数が増えてきてくれてようやく4人になった。かなり内輪感があふれるAdvent Calenderだと言う感じ。しかも内容もニッチなところを攻めている記事ばかり。そんなでも読んでくれる方には感謝しかない。

一応私の回では4回分(今回は3回目)の連載を考えているので、一通り目を通してくださると何となくストーリーが繋がると思う。

目標

Network上の拡散方程式を定義してその性質を見る

前提知識

僕が書いた記事を読んでいてくれればよい。特にLaplacianの方は必須。

• Scale Free NetworkをMatlabで実装する
• NetworkのLaplacianについて

Network上の拡散方程式の定義

さて、前回はNetworkにはLaplacian行列なるものが定義できることを学び、Network上の拡散現象を表すのに使えるんだ、という話をした。今回はその拡散方程式を定義していこうと思う。

ではNetwork上の拡散方程式をモデリングしていこう。ここでいうモデリングとは数理モデリングのことを意味し、現象を何らかの仮定を起きながら数式で表していく行為のことを言う。

ではモデリングする上での仮定を列挙する。ここでは時間tに対して、次の時間t + dtにて何が起きるかを見ていく。

1. Networkの各Node上に「何か」があり、Edgeを介して隣のNodeへとその「何か」が拡散する。
2. 拡散の仕方はFickの法則(リンクはWikipedia)に従う。すなわち、「流速は濃度勾配に比例する。」
3. 各Nodeについて(時刻t + dtでの量) = (時刻tでの量) + (dtで増えた量)

以上の仮定を元に、簡単な例を使って考えていく。まずは仮定1と2を見よう。
下の図のように、4つのNodeからなるNetworkの各Nodeの上に四角い何かがある状況を考える。これが仮定2に従って拡散するとはどういうことかを考える。

時刻ｔの状態

「流速」とはざっくり言うとものが流れる度合い。
「濃度勾配」とは隣との濃度(ここでは四角の数)の差。
これら2つが比例するというのだから、今得体の知れない「流速」を「濃度勾配」のなんとか倍で表すことができる。簡単のため比例定数を1とすれば、両者は一致するので、dtの間に「何か」がどれだけ流れるかは下図のようになる。

Node1に着目してやれば、dtの間にはNode2から4つ、Node4からは3つ受け取ることになるので、時刻tでのNode1での量を$X_1(t)$とすれば、dtが経つ時のNode1の量の増加量というのは$$((X_2(t) - X_1(t)) + (X_4(t) - X_1(t))$$と表すことができる。同じ感じで他のNodeがどれだけ増えるかということも考えればよい。

では、最後に仮定3について考えていこう。Node1について
(t + dtでの量) = (tでの量) + (dtでの増加量)
としてモデリングするのだったから、$$X_1(t + dt) = X_1(t) + D dt((X_2(t) - X_1(t)) + (X_4(t) - X_1(t))$$と表せる。ここで、比例定数を1としていたのを一般の$D$(所謂拡散係数)に戻し、dtも比例定数からひねり出した。では、式変形をしよう。右辺の$X_1(t)$を左辺に持ってきて、dtで全体を割って、右辺のかっこの中身を整理すれば、
$$\frac{X_1(t + dt) - X_1(t)}{dt} = D (-2X_1(t) + X_2(t) + X_4(t))$$
となる。左辺に注目すると、なんだか微分の定義式に見えてくるので、$dt \to 0$とすれば
$$\frac{dX_1}{dt} = D (-2X_1 + X_2 + X_4)$$
といった常微分方程式ができあがる。これがNode1に関する方程式である。同様にして、それぞれのNodeについても考えてあげれば、

\begin{eqnarray*}
  \begin{cases}
    \frac{dX_1}{dt} = D (-2X_1 + X_2 + X_4) & \\
    \frac{dX_2}{dt} = D (X_1 - 3X_2 + X_3 + X_4) & \\
    \frac{dX_3}{dt} = D (X_2 - X_3) & \\
    \frac{dX_4}{dt} = D (X_1 + X_2 - 2 X_4) &
  \end{cases}
\end{eqnarray*}

これを行列をつかってまとめてあげると、

\frac{d}{dt}\left(
\begin{matrix}
 X_1 \\ X_2 \\ X_3 \\ X_4
\end{matrix}\right) = D
\left(
\begin{matrix}
 -2 & 1 & 0 & 1 \\
 1 & -3 & 1 & 1 \\
 0 & 1 & -1 & 0 \\
 1 & 1 & 0 & -2
\end{matrix}\right)
\left(
\begin{matrix}
 X_1 \\ X_2 \\ X_3 \\ X_4
\end{matrix}\right)

となる。これがNetwork上の拡散方程式である。右辺の行列をじーっと見ると、何とこの行列はNetworkのLaplacian行列になっているのである。

つまり、上のことを一般化してあげれば、あるNetworkが与えられたとき、そのLaplacian行列を$L$、拡散係数を$D$、各Nodeの量を上から並べてあげたベクトルを$\boldsymbol{X}$とすれば、$$\frac{d\boldsymbol{X}}{dt} = DL\boldsymbol{X}$$とすっきりした形でNetwork上の拡散方程式を表すことができる。

「普通の」拡散方程式$$\frac{\partial u}{\partial t} = D\Delta u$$と見比べれば何故Laplacian行列が「Laplacian」と呼ばれるのか理解していただけるであろう。

Network上の拡散方程式の性質

以下かなり数学を使う。結果やSimulationだけを見たい人は次節へ進むことをおすすめする。

一様化

普通の拡散方程式は時間発展させてみると、ちゃんと空間に対して一様になってくれる。(少し騙している。境界条件をある条件にすればこれは正解。)水を張ったバケツにインクをぽちゃんとしたら、時間とともに一様に拡散していく様子を思い浮かべてみればわかると思う。これと同じように、Nodeに対して一様になってくれるような働きがきちんとNetwork上の拡散方程式にもあるかどうかを考察しよう。

何を確認すればよいかというと、拡散方程式$$\frac{d\boldsymbol{X}}{dt} = DL\boldsymbol{X}$$が時間発展の後に$\boldsymbol{X}$を一様化するかどうかを確認すればよい。一応簡単のため、Networkは強連結であることを仮定する。幸い、前回の通り、$L$の固有値は全て0以下であることから、このことを調べるのはそこまで難しくない。ちなみに強連結性がどこで効いてくるかというと、0固有値が1つしか出ないというところで効いてくる。

まずは、$L$の固有ベクトルにはどんなものがあるかを確かめる。0固有値に対する固有ベクトルは$$\boldsymbol{v}_1 = (1, \cdots, 1)^T$$である。ここで$\boldsymbol{v}$の添字はLaplacianの固有値を大きい順に並べた時のインデックスである。

では、このことを踏まえた上で、拡散方程式の解を見てみると、常微分方程式の解を求めれば良いから、$\lambda^{[i]}$を$i$番目の固有値として、
$$\boldsymbol{X}(t) = \sum_{i = 1}^N C_ie^{D\lambda^{[i]}t}\boldsymbol{v}_i$$
と表せる。$C_i$は積分定数である。

この解が実際に$t \to \infty$としてときどうなるかは、$\lambda^{[i]}$が$i = 2, 3, \cdots, N$に対して0より小さくなるから、上の和の$i = 1$のときすなわち、0固有値の時のみが残る。つまり、
$$\lim_{t\to\infty}\boldsymbol{X}(t) = C_1\boldsymbol{v}_1 = C_1(1, \cdots, 1)^T$$
となることがわかる。どんな初期値から始めても時間が経てば、全てのNodeの値が同じになるということを意味しているので、目的は果たせた。

(一応、強連結性を仮定しない場合は、連結成分が複数出る場合を考えるのだが、この時は、連結部分で一様化されることが計算できる。つながっていない部分に対してはどう頑張っても拡散はしないことを想像してもらえれば良い。)

一様化のスピード

これから話すことは普通の拡散方程式では見られないことである。まず、Scale Free NetworkのLaplacian行列の性質として、固有ベクトルの成分が局在していることを思い出してもらいたい。特に小さい固有値に関する固有ベクトルを見ると、これは次数の多いNodeに対応する成分が局在して存在しているという性質があった。固有値が小さければ小さいほど、上の拡散方程式の解の指数関数部分を早く0へと収束させる。つまり、次数の多いNodeほど、早く一様な状態になるということである。

Simulation

実際にNetwork上の拡散方程式を実装してみよう。例の如くmatlabで行う。Networkは以前実装したScale Free Networkにて行う。

まずは400個のNodeに適当な初期状態をばらまく。真ん中にあるほど次数が多くなるような描画方法を取っている。

時間経過。次数の多いところから緑(大体5)になっていく。

ほぼ定常状態。全てのNodeが大体5になった。

Simulationは4次のRunge-Kutta法を使用した。以下にコードを記載する。scalefreenetwork関数はここに載っている。

networkdiffusion.m

startSize = 2;
m = 11;
N = 400;
D = 0.06;
dt = 0.05;
randWidth = 10;
epsilon = 0.01;

graph = scalefreenetwork(startSize, m, N);
ave = mean(degree(graph));
x = randWidth * rand(N, 1);

colormap(jet);

L = - laplacian(graph);
e = eig(L);

for t = 0 : 1000
    prevX = x;
    x = nextXRK(x, dt, L, D);

    plot(graph, 'LineWidth', 0.2, 'EdgeAlpha', 0.1, 'MarkerSize', 10, ...
    'NodeCData', x, 'Layout', 'force');
    xticks([]);
    yticks([]);
    
    caxis([0, randWidth]);
    time = strcat("t = ", string(t * dt), "       1");
    time = extractBefore(time, 12);
    xlabel(time, 'FontName', 'FixedWidth', 'FontWeight', 'Bold');
    colorbar    

    pause(0.00001);
end

% this calculates next x with Runge-Kutta
function nextX = nextXRK(x, dt, L, D)
    k1 = D * L * x;
    k2 = D * L * (x + dt / 2 * k1);
    k3 = D * L * (x + dt / 2 * k2);
    k4 = D * L * (x + dt * k3);
    
    nextX = x + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
end

結論

Network上の拡散方程式を定義して、一様化と一様化のスピードについて見た。