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

前回はBurgers方程式をRで計算したが、Burgers方程式の非線形項 $u u_{x}$ を線形項 $cu_{x}$ にした式は"移流拡散方程式"と呼ばれる。今回は移流拡散方程式をRで計算し、Burgers方程式とはまた違った解の挙動となることを視覚的に確認しようと思う。

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

移流拡散方程式とは

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

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

この左辺第2項の移流項は、Burgers方程式とは異なり、波の速さが波高によらず同じであることを表す。また、右辺の拡散項は、時間の経過と共に波が広がっていくことを表す。なお、 $c$ 、 $\nu (\gt 0)$ は定数。

Burgers方程式では、時間の経過と共に"衝撃波"と呼ばれる突き立った波が生じ、その後、次第に衝撃破が減衰していくことを見たが、移流拡散方程式では、波は一定の速さで進むため、波が一定の速さを保ったまま、次第に減衰していくことが想像される。

移流拡散方程式についても、関数変換により拡散方程式に帰着させることや、いくつかの方法で解析解を構成することができるが、この記事では立ち入らない。

移流拡散方程式の差分方程式

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

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

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

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

を得る。ここで、 $A= \nu \frac{ \Delta t}{\Delta x^2}$ 、 $\tilde{B}=c \frac{ \Delta t}{2 \Delta x}$ である。（Burgers方程式とは $\tilde{B}$ が異なるのみ。）

計算結果

$\nu=0.02$ 、 $c=1.00$ とし、初期条件 $(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

計算結果の考察

コサインカーブから始まるが、波高によらず流速は一定のため、一定の速さを保ったまま、前後に裾野を広げながら波高が減少していく様子が見て取れる。最後には水平面に帰着する。名前の通り、移流しながら拡散していく様子が確認できる。（なお、 $c=-1.00$ で計算すると左側に進む波になる。）

やはりBurgers方程式とは異なる動きである。

計算コード（Rコード）

こちらは計算に使ったコード。特に洗練されたものではないので参考までに。Burgers方程式のコードと係数および差分方程式の箇所を除いて全く同じ。

nu=0.02
c=1.00
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
B.tilda=c*dt/(2*dx)
for (n in 1:tt){
        result[1,n+1]=(A-B.tilda)*result[2,n]+(1-2*A)*result[1,n]+(A+B.tilda)*result[xx,n]
        for (m in 2:(xx-1)){
                result[m,n+1]=(A-B.tilda)*result[m+1,n]+(1-2*A)*result[m,n]+(A+B.tilda)*result[m-1,n]
        }
        result[xx,n+1]=(A-B.tilda)*result[1,n]+(1-2*A)*result[xx,n]+(A+B.tilda)*result[xx-1,n]
}

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

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