読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

deSolve パッケージで拡散方程式(熱伝導方程式)を数値的に解く

f:id:abrahamcow:20141124141446p:plain

溶質が均一に分布している材料 A から, 材料 B に溶質が供給されていく状況を考える.

一次元の拡散方程式は,
 \displaystyle \frac{\partial u}{\partial t}  = D \frac{\partial ^2 u}{\partial x^2}  \tag{1}

ここで, D は拡散係数と呼ばれる定数 (D > 0) である.

初期条件と境界条件を以下のように与える.

初期条件:
\displaystyle \left\{ \begin{array}{ll} u = 1 & (x \le 30) \\ u = 0 & (x > 30) \end{array} \right.

境界条件
\displaystyle u(x=0,t) = u(x=60,t)=0

この微分方程式を deSolve パッケージで解くには, (1) 式を

\displaystyle \left\{ \begin{array}{l} \partial u =D \frac{\partial u_1}{\partial x}    \\ u_1 = \frac{\partial u_2}{\partial x}  \end{array} \right.

と, 書き換える.

微分を差分で近似し,

\displaystyle \frac{du_1}{dt} = D \frac{u_{i,i+1 } -u_{i-1,i}}{\Delta x_i} ,

\displaystyle u_{i-1,i} = D \frac{u_{i} -u_{i-1}}{\Delta x_{i-1,i}} ,

この式をコードに与えてやる.

library("deSolve")
diffusion <- function(t, dens, parameters) {
  deltax <- c (0.5, rep(1, numboxes - 1), 0.5)
  #sum(deltax) 60
  du1 <- D * diff(c(0, dens, 0)) / deltax
  du2 <- diff(du1) / delx 
  # the return value
  list(du2)
} # end
D <- 0.3 #diffusion coefficient
delx <- 1 #thickness of boxes
numboxes <- 60
u <- rep(0, times = numboxes)
u[1:30] <- 1
state <- c(u = u)
times <-seq(0, 200, by = 1)
#print(system.time(
  out <- ode.1D(state, times, diffusion, parms = 0, nspec = 1, names="density")
#))

計算結果を下に図示する.

library(ggplot2)
library(tidyr)

colnames(out) <- c(colnames(out)[1],as.character(1:60))
out2 <- gather(as.data.frame(out),x, density, -time)
out2$x <-as.numeric(out2$x)

theme_set(theme_grey(20))
ggplot(out2,aes(x=x,y=density,group=time,colour=time)) +
  geom_line()
length(Distance)
sum(out2$time==0)

f:id:abrahamcow:20150219115637p:plain


参考文献:
http://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf

関連エントリ:
拡散方程式の手短な導出 - 廿TT
一次元の拡散方程式の解 - 廿TT