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

廿TT

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

一次元の拡散方程式の解

拡散方程式の解


拡散方程式の手短な導出 - 廿TT の続き.

f:id:abrahamcow:20141124141446p:plain

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

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

ここで,

  • c: 濃度
  • x: 位置
  • t: 時間
  • D: 拡散係数と呼ばれる定数 (D > 0)

である.

濃度 c は時間によらず一定(定常状態), すなわち \partial c / \partial t =0 とすると, (1) の解は,
 \displaystyle c = Ax+B \tag{2}
A, B積分定数)である。

全区間での濃度は一定
 \displaystyle \int^{\infty}_{-\infty} c(x,t) dx = M
M:単位面積あたりの溶質の全量)とすると,
 \displaystyle c(x,t)= \frac{M}{2 \sqrt{\pi D t}} \exp \left(-\frac{x^2}{4 D t} \right) \tag{2}
はこれを満たす.

なぜならば,  \int^{\infty}_{- \infty} e ^{-x^2} dx = \sqrt{\pi} はすでにご承知おきのこととして(ガウス積分 - Wikipedia),
 \displaystyle \int^{\infty}_{- \infty} e ^{-a x^2} dx = \frac{1}{\sqrt{a}} \int^{\infty}_{- \infty} e^{(-\sqrt{a}x)^2} d (\sqrt{a}{x}) = \sqrt{\frac{\pi}{a}}
だからである.

f:id:abrahamcow:20141124141852p:plain

さて, 初期条件が  x \le 0 c=c_0,  x>0 c=0 の場合を考える.

時刻 t において微小区間 dy の溶質の量は c_0 dy である.
よって, 時刻 t において距離 y だけ離れた, 位置 x にある点の濃度に寄与する量は, (2) 式において  M=c_0 dy と置き換えて,
 \displaystyle \frac{c_0 dy}{2 \sqrt{\pi D t}} \exp \left(-\frac{x^2}{4 D t} \right)
で与えられる.

 x < 0 の区間に  c_0 dy の溶質があるので, 積分範囲は  y >x である.

 \displaystyle c(x,t) =\frac{c_0}{2 \sqrt{\pi D t}} \int^{\infty}_{x} \exp \left(-\frac{x^2}{4 D t} \right) dy \\
\displaystyle \frac{c_0}{2 \sqrt{\pi}} \int^{\infty}_{x/ \sqrt{4 D t}} \exp^{-u^2} du

これは補誤差関数 erfc を用いて,
 \displaystyle c(x,t) = \frac{c_0}{2}\rm{erfc}\left(-\frac{x}{2 \sqrt{D t}} \right)
と表せる.

R による計算例

初期濃度を 1, 拡散係数を 10 と置いた場合, 拡散は下図のように進行する.

f:id:abrahamcow:20141124141423p:plain

erfc <- function(x) 2 * pnorm(x * sqrt(2), lower=FALSE)
c_func <- function(x, t=t1, c_0=1, D=10){
  (c_0/2) * erfc(x/(2 * sqrt(D*t)))
}

x_seq <- seq(-5,5,by=0.01)

mat1 <- matrix(,nrow=1000,ncol=length(x_seq))

length(x_seq)
t1=0
for(i in 1:1000){
  mat1[i,] <- sapply(x_seq,c_func)
  t1 = t1 +1
}
#nihongo()

plot(x_seq,mat1[1,], type="l",col="red3",lwd=2, xlab="x:位置", ylab="c:濃度")
lines(x_seq,mat1[2,], type="l",col="red3",lty=2)
lines(x_seq,mat1[1000,], type="l",col="red3",lty=3)
legend("topright", legend=c("t=1","t=2","t=1000"), lty=1:3, col="red3")

GIF アニメも作ったんだが, ループさせるやり方がわからなかった.
新しいタブで画像を開きなおしてリロードして再生してください.
f:id:abrahamcow:20141124142145g:plain

library(animation)
saveGIF({
  par(family = "HiraKakuPro-W3")
  for(i in 1:100){
    plot(x_seq,mat1[i,], type="l",col="red3",lwd=2, xlab="x:位置", ylab="c:濃度", ylim=c(0,1))
  }
},interval=0.5)