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

廿TT

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

拡散方程式, 差分法, 陽解法, Rによる練習.

目的

拡散方程式の差分解法として最も簡単と思われる陽解法をとりあえず動かしてみる.

※本文中のコード、あってるか自信ない.

拡散方程式の差分表式

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

ここで,

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

である.

今回は a = 1 とさせてもらう.

u(x,t) を 単に u と省略して書くことにする.

 \displaystyle \frac{\partial u}{\partial t} = \frac{\partial ^2 u}{\partial x^2} \tag{1}

この方程式を.
初期条件:u(x=0) = φ(x), (0 < x < 1 )
境界条件u(0,t) = 0, u(1, t) = 0, ( t > 0 )
のもとで解く.

今回初期値は,
 \phi(x) =\left\{
\begin{array}{ll}
x & 0\le x \le1/2\\
1-x & 1/2 \le x \le 1
\end{array}
\right.

という折れ線にする.

(1) 式の微分を差分で近似し,

\displaystyle \frac{u^{(n+1)}_{j} - u^{(n)}_{j}}{\Delta t} =
\displaystyle \frac{( u^{(n)}_{j+1} - u^{(n) }_{j} ) - (u^{(n) }_{j} - u^{(n) }_{j-1})}{(\Delta x)^2} \\
\displaystyle =\frac{( u^{(n)}_{j+1} - 2u^{(n) }_{j}   + u^{(n) }_{j-1})}{(\Delta x)^2}

という方程式を得る.

これを u^{(n)}_{j+1} について解き, 初期条件と境界条件も併せて書くと, 以下の更新式を得る.

\displaystyle u^{(n+1)}_{j} =u^{(n)}_{j} + \frac{\Delta t}{(\Delta x)^2} (u^{(n)}_{j+1} -2 u^{(n)}_{j}+ u^{(n)}_{j+1} )
(j=1,2, \ldots, J-1; \, n= 0,1, 2, \ldots),

u^{(0)}_{j} = \phi (j \Delta x) \quad (j=1,2, \ldots, J-1) ,

u^{(n)}_{0} =0, ~ u^{(n)}_{J} =0 \quad (n=0,1,2, \ldots) .

こうやって微分方程式を数値的に解く方法のことを陽解法という.

Rによる練習

f:id:abrahamcow:20150201172224p:plain

 \Delta x = 1/8, \Delta t = 1/160 とした.

上から順に  j = 0,1,2 各時刻の u を表している.

phi <- function(x){
  ifelse(x < 1/2,x,1-x)
}

Deltax <-1/8
jDeltax <- seq(0,1,by=Deltax)
J <- length(jDeltax) 

Deltat <-1/160

u_old <- phi(jDeltax)
plot(jDeltax,phi(u_old),type="b",lwd=1.5, lty=1, xlab="x", ylab="u")
tmp <- numeric(J)
#ここで for を使うのは R らしくない気がする...
for(i in 1:2){
  for(j in 2:(J-2)){
      tmp[j] <- 
      u_old[j+1] -  2*u_old[j] + u_old[j-1] 
  }
  u_new <- u_old + (Deltat/(Deltax^2)) * (tmp)
  points(jDeltax,u_new,type="b",lwd=1.5, lty=(i+1))
  u_old <- u_new
}

参考文献

偏微分方程式の差分解法 (東京大学基礎工学双書)

偏微分方程式の差分解法 (東京大学基礎工学双書)