Rで偏微分方程式 - 餡子付゛録゛

速度的に実用性が低い感じもするのですが、Rで偏微分方程式の数値解を計算することもできます。微分方程式一般に言えることですが、せっせと計算して閉じた解を求めるのは苦労するので、数値解は実用的には重要です。
deSolveパッケージを使えば三元連立式ぐらいまでは容易で、使い方はサンプル・コードを見ればよいのですが、恐らく最も見慣れている偏微分方程式である波動方程式のサンプルを書いてみました。
deSolveの旧バージョンであるReacTranパッケージのリファレンスに例としては上がっているわけですが、パッケージの構成が変わってそのままでは動かないので、もしかしたら参考になるかも知れません。

1. 計算の概要

まずは波動方程式を眺めてみましょう。波動関数は $u(t,x)$ で、 $x$ を位置、 $t$ を時間とします。 $c$ は波動の速度。
$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$
偏微分方程式のままだと解けないので、常微分方程式に変形します。
$\frac{\partial u_1}{\partial t} = u_2$
$\frac{\partial u_2}{\partial t} = c^2 \frac{\partial^2 u_1}{\partial x^2}$
初期パラメーターを以下のようにします。力学的振動のときと異なって、ベクトルで大量の値がセットされることに注意してください。
$u_1(0,x)=exp{-0.05x^2}$
$u_2(0,x)=u(t,\infty)=u(t,-\infty)=0$
$c=1$
この式からそのままコーディングできる人は、かなり数値演算に慣れていると思います。

2. 実際の演算コード

以下が実際のコードになります。ポイントは恐らく、二階微分を行うところになると思います。数値微分のちょっとした知識が要ります。deSolveパッケージの癖で、複数の変数のベクトルが、一つのベクトルにまとめられているのも、注意しないと混乱するかも知れません。

# deSolveパッケージを読み出す
library("deSolve")
# 位置xの刻み幅は0.2
dx <- 0.2
# xの最小値と最大値
xlim <- c(-70, 70)
# -40から40まで位置xのベクトルを得る
x <- seq(xlim[1], xlim[2], dx)
# xの数を保存しておく
N <- length(x)
# 波の高さの初期値
uini <- exp(-0.05 * x^2)
# 波の高さの変化分の初期値
vini <- rep(0, N)
# パラメーターは同じベクターにしておく
yini <- c(uini, vini)
# 時間tのベクトルを得る
times <- seq (from = 0, to = 50, by = 1)
# 常微分方程式の形式で、波動方程式を表現する
wave <- function(t, y, parms, dx, N) {
# ベクターを二つに分ける
u1 <- y[1:N]
u2 <- y[(N+1):(2*N)]
# 波の高さの変化を計算
du1 <- u2
# 波の高さの変化分の変化の計算
# 二階微分を行う。両端の変化率は境界条件によりゼロ。
du2 <- c(0, (diff(u1[-1],1) - diff(u1[-N],1))/dx^2, 0)
# 変化を戻り値で返す
return(list(c(du1, du2)))
}
# ode.1D関数は一次元の微分方程式を扱う。
# nspecは計算する変数の数で、ここではu1とu2で二つ。
# 計算手法はデフォルトだと収束しなくなるので、ode45
out <- ode.1D(yini, times, wave, parms = 0, nspec = 2, dx=dx, N=N, method="ode45", atol=1e-14, rtol=1e-14)

ちょっと時間がかかると思いますが、気長に待てば計算できるでしょう。

3. 演算結果をプロットしてみる

出てきた計算結果は以下のようにプロットしたり、中身を見ることができます。

# 波の変化をプロットしてみる
matplot.1D(out, which=1, subset = time %in% seq(0, 50, by = 10), type = "l", col = c("black", rep("darkgrey", 5)), lwd = 2, grid = x, xlim = xlim, main="u(x, t)")
legend("topright", lty = 1:6, lwd = 2, col = c("black", rep("darkgrey", 5)), legend = paste("t = ",seq(0, 50, by = 10)))

# 等高線プロットをしてみる
image(out, xlab = "t", ylab = "x", main = "PDE", add.contour = FALSE)

左が $u(t, x)$ 、右が $du$ になります。色が波動の高さで、赤色が高く青色が低くなります。表題などを綺麗にする方法は調査中。

# t=30のときの大きさを波の形状を見てみる。ylim重要。
plot(out[30,][1:N], ylim=c(0, 1), xlab="x", ylab="u")
# 中央地点x=2*N/3の波の高さの変化を見る
plot(out[,2*N/3][1:length(times)], ylim=c(0, 1), xlab="t", ylab="u")

コンピューターを使っている気になれますね。