グラフ上の離散化された拡散方程式
準備
離散化された拡散方程式の状態空間モデル の続き的な内容.
1次元→2次元→3次元と拡張していってもいいのだけれど,ここではグラフ(ノードやエッジがある方のグラフ)上の拡散を考えてみる.
さて,1次元の拡散方程式グラフ上
を離散化し,
この差分方程式を次の画像のようにグラフで描いて考えてみる.
グラフの丸をノード,線をエッジと呼ぶことにする.
(1)式の右辺の
一方,
変数を行列・ベクトルでまとめて次の図のように表せる.
結果,時点
は,隣接関係を表す行列(Adjacency),エッジの数を表す行列(Degree)と単位行列(Unit)を使って,
のように表せることがわかった.
これをグラフ上の拡散の定義だと思って,隣接関係を次の図のような上下左右として数値計算してみよう.
Rによる数値例
といってもちょっと力尽きてしまって,いくつか図を描くだけである.
隣接関係を表す行列を作る関数は次のようにした.
make_adjmat_2d <-function(nr, nc){
nn = as.integer(nr*nc)
A <- Matrix::spMatrix(nrow = nn, ncol = nn)
edges = nr*(0:nc)
for(i in 1L:nn){
if(i-nr > 0){ #left bound
A[i,i-nr] <- 1
}
if( i %% nr != 1L ){ #top bound
A[i,i-1] <- 1
}
if((i+nr) <= nn){ #right bound
A[i,i+nr] <- 1
}
if( (i %% nr != 0L) & i+1L <= nn ){ #bottom bound
A[i,i+1] <- 1
}
}
return(A)
}
if 文を使って上下左右の端っこのところで隣接関係がないようにしている.
この行列をプロットするとこんなふうだ:
1期先の状態を作るのは次のようにした.
difit <- function(coef, u, nt){
U <- matrix(0, nt+1, length(u))
U[1,] <- u
for(i in 1:nt){
u2 <- coef%*%u
U[i+1,] <- as.vector(u2)
u <- u2
}
return(U)
}
これは nt
ステップ文 for 文で回すだけであとは上の式そのままという感じ.
係数行列(上の式の
makecoef_exp <- function(r, A){
diag(nrow = nrow(A)) + r*(A-diag(rowSums(A)))
}
これも上の式そのままという感じ.
この行列をプロットするとこんなふうだ:
初期状態で2箇所だけ濃度が高いところを作って24ステップ回した.
uini <- numeric(length = nrow(A))
uini[1] <- 1
uini[8] <- 1
U <- difit(beta, uini, 24)
print(rowSums(U))
次のように,全体での量が変化していないことをたしかめた.
> print(rowSums(U))
[1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
拡散していく様子のGIFアニメ:
いくつか切り口を変えた図も描いてみる.
とりあえず全部並べる:
横(x)方向に切り取って濃度の折れ線グラフ:
縦(y)方向に切り取って濃度の折れ線グラフ:
最後に,使用したRのコード全体はこちら:
Discussion