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

前回は移流拡散方程式をRで計算したが、移流拡散方程式の右辺を0としたもの、つまり移流項のみとした式は"移流方程式"と呼ばれる。今回は移流方程式をRで計算し、解の挙動を視覚的に確認しようと思う。（扱う方程式がだんだん易しくなっている）

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

移流方程式とは

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

$\dfrac{\partial u}{\partial t} + c\dfrac{\partial u}{\partial x} =0$

非常にシンプルな式であり、よっぽど特殊なケースを除けば、最も易しい偏微分方程式なのではないかと思われる。

さて、この解を波と考えることにすると、この解である波については、定数係数の移流項により、波高によらず波が一定の速さで進んでいくこと、および、拡散項が無いので、波の減衰は生じないことが想像される。つまり、初期波形を保ったまま同じ速さで進んでいくことが想像される。

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

移流方程式の差分方程式

続いて、数値計算用に移流方程式を差分方程式に書き換える。ここでも、時間微分に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} = 0$

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

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

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

計算結果

$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$ で計算すると左側に進む波になる。）

予想通りである。

計算コード（Rコード）

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

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))の下で計算
B.tilda=c*dt/(2*dx)
for (n in 1:tt){
        result[1,n+1]=result[1,n]-B.tilda*(result[2,n]-result[xx,n])
        for (m in 2:(xx-1)){
                result[m,n+1]=result[m,n]-B.tilda*(result[m+1,n]-result[m-1,n])
        }
        result[xx,n+1]=result[xx,n]-B.tilda*(result[1,n]-result[xx-1,n])
}

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

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