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

前回はKdV方程式をRで計算したので、今回はBurgers方程式をRで計算し、KdV方程式とはまた違った解の挙動となることを視覚的に確認しようと思う。なお、バーガースはオランダの物理学者(1895－1981)。

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

Burgers方程式とは

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

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

この左辺第2項は非線形項であり、波の速さが波高によって異なることを表す。KdV方程式にも見られた項であり、移流項と呼ばる。一方、右辺は $u$ の二階微分であり、三階微分だったKdV方程式とは異なる項である。右辺のこの項は拡散項と呼ばれ、時間の経過と共に波が広がっていくことを表す。なお、定数 $\nu \gt 0$ は動粘性係数と呼ばれる。

また、Burgers方程式は、 $u=-2\nu \frac{\partial \log \theta}{\partial x}$ と変換（この変換はCole-Hopf変換と呼ばれる）すると、 $\frac{\partial \theta}{\partial t}=\nu \frac{\partial^2 \theta}{\partial x^2}$ となって、拡散方程式（熱伝導方程式）に帰着される。そのため、ここからBurgers方程式の解析解を構成することができるが、この記事では立ち入らない。

Burgers方程式の差分方程式

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

$\dfrac{u_{m}^{n+1}-u_{m}^{n}}{\Delta t}+u_{m}^{n} \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} \\=\left( A - B \right)u_{m+1}^{n} + \left( 1- 2A \right)u_{m}^{n}+\left( A+B \right)u_{m-1}^{n}$

を得る。ここで、 $A= \nu \frac{ \Delta t}{\Delta x^2}$ 、 $B=u_{m}^{n} \frac{ \Delta t}{2 \Delta x}$ である。

計算結果

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

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

www.youtube.com

計算結果の考察

コサインカーブから始まるが、波高が高いほど流速が速いため、時間の経過と共に波が突き立つ（衝撃波と呼ばれる）。しかしながら、突き立った波は次第に減衰していき、ついには水平面に帰着する。

KdV方程式の解とは異なり、Burgers方程式の解はある意味常識的な挙動を示す。方程式上は $u$ の三階微分か二階微分かの差だが、記述される現象が大きく異なるのが興味深い*1。

計算コード（Rコード）

こちらは計算に使ったコード。特に洗練されたものではないので参考までに。

KdV方程式の計算コードと枠組みはほとんど同じで、計算結果を格納するResultという行列（x軸を行、t軸を列）について、時間0の初期値を入れるために最初の1列のみ（tt=1）を使用した点が異なる。（KdV方程式では初期値の格納に最初の2列を使用した。）

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

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

Burgers方程式についても、かつてCで計算したことがありましたが、今回Rで計算してみました。間違い等ありましたらお手数ですがコメント等いただけますと幸甚です。

*1:では右辺が $u$ の一階微分ならどうなるかと思うかもしれない。 $u_{t}+uu_{x}=\nu u_{x}$ の式は、 $u-\nu$ をあらためて $u$ とすれば、 $u_{t}+uu_{x}=0$ となるので、Burgers方程式で $\nu=0$ とすることと本質的に同じである。そしてこの方程式でも当初は衝撃波が生じるが、減衰しないため、多値となって計算不能になる。