拡散方程式のRによる計算 - 吾未知足唯修身爾

前回は移流拡散方程式から拡散項を外した「移流方程式」をRで計算したが、移流拡散方程式から移流項を外した式、つまり拡散項のみとした式は、"拡散方程式"と呼ばれる。なお、この式は、"熱伝導方程式"としても知られるとても有名な式である。今回は、この拡散方程式をRで計算し、解の挙動を視覚的に確認しようと思う。

拡散方程式とは
拡散方程式の差分方程式
計算結果
計算結果の考察
計算コード（Rコード）

拡散方程式とは

拡散方程式とは、次の形の線形偏微分方程式である。

$\dfrac{\partial u}{\partial t} =\nu \dfrac{\partial^{2} u}{\partial x^{2}}$

右辺の拡散項は、時間の経過と共に初期値が広がっていくことを表す。なお、 $\nu (\gt 0)$ は定数。

解（例えばこれを波と考えるとすると）は、移流項が無いので時間が経っても一切進まず、一方で、拡散項により時間と共に波高が減衰していくだろうことが想像される。つまり、初期波形が減衰しながら水平面に帰すことが想像される。

なお、例によって、解析解についてはこの記事では立ち入らない。

拡散方程式の差分方程式

続いて、数値計算用に拡散方程式を差分方程式に書き換える。ここでも、時間微分に1次の前進差分、空間微分に2次の中心差分を用いることとすると、

$\dfrac{u_{m}^{n+1}-u_{m}^{n}}{\Delta t}= \nu \dfrac{u_{m+1}^{n}-2u_{m}^{n}+u_{m-1}^{n}}{\Delta x^2}$

となるので、これを整理すると、

$u_{m}^{n+1} =A u_{m+1}^{n} + ( 1- 2A )u_{m}^{n}+ A u_{m-1}^{n}$

を得る。ここで、 $A= \nu \frac{ \Delta t}{\Delta x^2}$ である。（移流拡散方程式で $\tilde{B}=0$ としたものと同じ。）

計算結果

$\nu=0.03$ とし、初期条件 $(1.0+\cos \pi x)$ および周期境界条件 $u(0,t)=u(2,t)$ のもと上記の差分方程式に従って計算する。ここで、 $dx=0.01$ 、 $dt=5.0\times10^{-4}$ とする。（これまでのKdV方程式、Burgers方程式、移流拡散方程式、移流方程式と全く同じ初期条件・境界条件、メッシュサイズ）

解の挙動は以下の通り。今回も一時停止や巻戻しが出来るようyoutube経由。

www.youtube.com

計算結果の考察

コサインカーブから始まるが、流速はゼロのため、初期波形が右（ないし左）に進むことはない。一方で、時間の経過とともに波高は減衰し、最終的には水平面に帰す様子が見て取れる。

計算コード（Rコード）

こちらは計算に使ったコード。特に洗練されたものではないので参考までに。移流拡散方程式のコードで $\tilde{B}=0$ としたものに相当する。

nu=0.03
dx=0.01
dt=5*10^-4

#計算範囲の指定(x=0~2,t=0~10)
xx=2/dx
tt=10/dt
result=matrix(0, nrow=xx, ncol=tt+1)

#初期条件の指定(1+cosπx)
for (m in 1:xx){
        result[m,1]=1+cos(pi*dx*(m-1))
}

#周期境界条件(U(0,t)=U(2,t))の下で計算
A=nu*dt/dx^2
for (n in 1:tt){
        result[1,n+1]=A*result[2,n]+(1-2*A)*result[1,n]+A*result[xx,n]
        for (m in 2:(xx-1)){
                result[m,n+1]=A*result[m+1,n]+(1-2*A)*result[m,n]+A*result[m-1,n]
        }
        result[xx,n+1]=A*result[1,n]+(1-2*A)*result[xx,n]+A*result[xx-1,n]
}

描画の方は省略。なお、面グラフの色は"orange"としている。

間違い等ありましたらお手数ですがコメント等いただけますと幸甚です。