数理解析研究所講究録
第 1800 巻 2012 年 226-235226

非線形波動方程式の安定かつ高精度な数値解法について

大阪府立大学工学研究科

北川 真帆 (Maho Kitagawa),

村上 洋一 (Youichi Murakami)

Graduate School of Engineering
Osaka Prefecture University
1 はじめにこの研究では，次のような移流型の非線形項を持つ 1 次元の偏微分方程式 (バーガース方程式
や $KdV$ 方程式等) の時間発展の数値解法について考察する．
$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+L(\frac{\partial}{\partial x},\frac{\partial^{2}}{\partial x^{2}},\frac{\partial^{3}}{\partial x^{3}},$

ここで，$L( \frac{\partial}{\partial x},\frac{\partial^{2}}{\partial x^{2}},\frac{\partial^{3}}{\partial x^{3}},$$\ldots)$(1)$\ldots)u=0$ ,は，微分の高階微分を含んだ線形演算子を示す．
$x$スペクトル法を用いて，時間発展方程式に適切な時間に関する差分法 (時間スキーム) を見出すことを目的とする．陽解法で数値安定になり精度も柔軟に選択できるものに着目する．

まず，取り扱う移流型の非線形項を持つ時間発展方程式の性質について簡単に述べ，次に数値計算法について説明する．特に，今まで取り扱われてきた時間差分法やこの研究で提案する方法を系統的にまとめる．バーガース方程式 (エネルギー非保存),$KdV$ 方程式(エネルギー保存)等について種々の時間差分法による数値安定条件の結果を与え，比較する．また，時間差分法による精度の違いを $KdV$ 方程式の保存量によって検討する．最後に，この研究の結論と今後の課題についてまとめる．

2 時間発展偏微分方程式
バーガース方程式$\frac{\partial u}{\partial t}=-\frac{1}{2}\frac{\partial u^{2}}{\partial x}+v\frac{\partial^{2}u}{\partial x^{2}}$,(2)は NavierStokes 方程式の簡単な 1 次元モデルとして Burgers[1] により導入された．ここで，$v$は正の定数である．非線形性と散逸性を持つ最も簡単な非線形拡散方程式である．
$KdV$ 方程式 [2],$\frac{\partial u}{\partial t}=-\frac{1}{2}\frac{\partial u^{2}}{\partial x}-\mu\frac{\partial^{3}u}{\partial x^{3}}$(3)は長波長の浅水波を記述する方程式として，非粘性流体の渦なしの方程式と境界条件から導かれた．ここで， は定数である．この方程式は無限の保存量を持つ保存系である．簡単なものを以$\mu$

下に 2 つ与える [3].

227(4)$M[u(t)]:= \int_{0}^{2\pi}u^{2}(t,x)dx$ ,$E[u(t)]:= \frac{1}{2}\int_{0}^{2\pi}(\mu(\partial_{\chi}u(t,x))^{2}-\frac{1}{3}u^{3}(t,x))dx$.

(5)$KdV-$ バーガース方程式 [4] とは次のような形で表される．
$\frac{\partial u}{\partial t}=-\frac{1}{2}\frac{\partial u^{2}}{\partial x}+v\frac{\partial^{2}u}{\partial x^{2}}-\mu\frac{\partial^{3}u}{\partial x^{3}}$

ここで，
$v,$

$\mu$

(6)

,は定数である．一般に分散性と散逸性を持っ非線形系においては，弱非線形を示す

小さな振幅の長波長の波を考えたときにこれら 2 つの効果が同程度のオーダーならば $KdV-$ バーガース方程式が導かれる．

5th $KdV$ 方程式は次式で与えられる．
$\frac{\partial u}{\partial t}=-\frac{1}{2}\frac{\partial u^{2}}{\partial x}+\gamma\frac{\partial^{5}u}{\partial x^{5}}$(7),

ここで，は定数である．一般の分散性を持つ非線形系においては，弱非線形を示す小さな振幅の$\gamma$$KdV$ 方程式が導かれるが多いが，まれに 3 階微分の係数がほとんど長波長の波を考えたとき，$0$になることがある．このような場合は高階の微分項である 5 階微分の分散項を考慮する必要がある．

3 数値スキーム
空間分割にはスペクトル法を用い，時間発展には 4 次のルンゲークッタ法，ジェイムソンの方法，積分因子法を用いる．また，2/3 則によってアライアジングエラーを除去する．それぞれの方法について詳しく説明する．

31 スペクトル法
スペクトル法とは，関数空間で直交関数展開して方程式を近似する方法である．この方法は周期境界条件の時以外では使用できないが，すべての成分を使って微分を近似しているため差分法に比べ精度がよい．

324 次のルンゲークッタ法
ここでは微分方程式 (8) に基づいて 4 次のルンゲークッタ法について説明する．
$\frac{dy}{dt}=f(t,y)$ある時亥$it_{i}\ovalbox{\tt\small REJECT}$.

(8)

こおける $y$ の値$y_{i}$がわかっているとき，$t_{i+l}=t_{j}+\Delta t$

における$y$の値$y_{i+1}$を次のように近似する．
$y_{n+1}=y_{n}+ \frac{\Delta t}{6}(k_{1}+2k_{2}+2k_{3}+k_{4})+O(\Delta t^{5})$,(9)ただし、$k_{1}=f(t_{n},y_{n})$

,

(10)

228
$k_{2}=f(t_{n}+ \frac{\Delta t}{2},y_{n}+\frac{\Delta t}{2}k_{1})$

,

(11)

$k_{3}=f(t_{n}+ \frac{\Delta t}{2},y_{n}+\frac{\Delta t}{2}k_{2})$

,

(12)

$k_{4}=f(t_{n}+\Delta t,y_{n}+\Delta tk_{3})$

.

(13)

5 次以上になるとルンゲークッタ法の式は非常に複雑になる．また，5 次のルンゲークッタ法は6 段必要となり，あまり使われない．

33 ジェイムソンの方法
ここでは微分方程式 (14) に基づいてジェイムソンの方法 [5] について説明する．ジェイムソンの方法はルンゲークッタ法の特別な方法である．
$\frac{dy}{dt}=f(y)$

.

(14)

ただし、 $f$ が にあらわに依存している関数の場合この方法は適用できない．
$t$

まず、 1 次のオーダーで考える。 式 (14) を差分化すると，$y_{n+1}=y_{n}+\Delta tf(y_{n})$$=y_{n}+ \Delta t\frac{dy_{n}}{dt}+0(\Delta t^{2})$(15)が得られる．次に，2 次のオーダーは，$y_{(1)}=y_{n}+ \frac{\Delta t}{2}f(y_{n})$$=y_{n}+ \frac{\Delta t}{2}\frac{dy_{n}}{dt}$

(16)

,

$y_{n+1}=y_{n}+\Delta tf(y_{(1)})$$=y_{n}+ \Delta tf(y_{n}+\frac{\Delta t}{2}\frac{dy_{n}}{dt})$$=y_{n}+ \Delta t\frac{d}{dt}(y_{n}+\frac{\Delta t}{2}\frac{dy_{\mathfrak{n}}}{dt})$$=y_{n}+ \Delta t\frac{dy_{n}}{dt}+\frac{\Delta t^{2}d^{2}y_{n}}{2dt^{2}}+0(\Delta t^{3})$(17)となる．さらに，ジェイムソンの方法は N 次のオーダーに一般化することが可能で，以下のように書くことができる．
$y_{(1)}=y_{n}+ \frac{\Delta t}{N}f(y_{n})$

(18)

,

$y_{(2)}=y_{n}+ \frac{\Delta t}{N-1}f(y_{(1)})$

,

(19)

229$y_{(N-1)}=y_{n}+ \frac{\Delta t}{2}f(y_{(N-2)})$

(20)

,

$y_{n+1}=y_{n}+\Delta tf(y_{(N-1)})+0(\Delta t^{N+1})$

.

(21)

ジエイムソンの方法は途中の値を保存する必要がないので，配列 (メモリー) が少なくて済む．
また，通常のルンゲークッタ法のように次数により段数が増えるなど複雑になることもない．
34 積分因子法
ここでは微分方程式 (22) に基づいて積分因子法 [6] について説明する．

.

$\frac{dy}{dt}=f(y)-ay$(22)

式 (22) を次のように式変形する．ここで，$f(y)$ は$\frac{dy}{dt}+ay=f(y)$$y$の非線形関数である．
(23)

,

(24)

$( \frac{dy}{dt}+ay)e^{at}=f(y)e^{\alpha t}$

すると，以下のようにまとめることができる．
$\frac{d}{dt}(e^{\alpha t}y)=e^{\alpha t}f(y)$.

(25)

この形の元で時間微分を差分化するのが積分因子法である．例えば時間について前進オイラー法
を用いると式 (26) のようになる．
$y_{n+1}=e^{-\alpha\Delta t}(y_{n}+\Delta tf(y_{n}))$.

積分因子法は線形項を解析的に解いているため，線形項のみなら$\Delta t$(26)$\Delta t$による打切り誤差が生じず，を大きくしたことによる数値不安定も生じない．この研究では積分因子法を効果的に用いて数値安定な数値計算法を見出すことを目的とする．

35 ルンゲークッタ法と積分因子法を組み合わせる方法
ここでは微分方程式 (27) に基づいて 4 次のルンゲークッタ法と積分因子法を組み合わせる方法にっいて説明する．
$\frac{dy}{dt}=f(y)-\alpha y$

式 (27) に積分因子$e^{\alpha t}$.

(27)をかけ，以下のように式変形する．

$\frac{d}{dt}(e^{\alpha t}y)=e^{at}f(y)$.

(28)

式 (28) に 4 次のルンゲークッタ法を用いると式 (29) のようになる．
,$y_{n+1}=e^{-\alpha\Delta t}(y_{n}+ \frac{\Delta t}{6}m_{1})+e^{-\frac{a}{2}\Delta t}\frac{\Delta t}{3}(m_{2}+m_{3})+\frac{\Delta t}{6}m_{4}+0(\Delta t^{5})$

(30)

230
ただし，$m_{1}=f(y_{n})$

(31)

,

$m_{2}=f(exp(- \frac{\alpha}{2}\Delta t)(y_{n}+\frac{\Delta t}{2}m_{1}))$$m_{3}=f(exp(- \frac{\alpha}{2}\Delta t)y_{n}+\frac{\Delta t}{2}m_{2})$

(32)

,

(33)

,

$m_{4}=f(exp(- \alpha\Delta t)y_{n}+\Delta texp(-\frac{a}{2}\Delta t)m_{3})$

.

(34)

35 線形項と非線形項を段階的に解く方法
ここでは線形項 -僻と非線形項 $f(y)$ を持つ微分方程式 (35) に基づいて線形項と非線形項を段階的に解く方法について説明する．
$\frac{dy}{dt}=f(y)-\alpha y$.

(35)

まず，積分因子法で線形項のみを時間発展させ，$\frac{d\tilde{y}}{dt}=-\alpha y$

(36)

,

$\tilde{y}=e^{-\alpha\Delta t}y_{n}$(37),この項を用いて非線形項を解く．その際，適当な次数のルンゲークッタ法やジェイムソンの方法を用いる．例えば，前進オイラー法を用いると次のようになる．
$\frac{dy}{dt}=f(\tilde{y})$(38),$y_{n+1}=\tilde{y}+\Delta tf(\tilde{y})$

.

(39)

4 結果
各方程式についてそれぞれ数値安定性を調べた．まず空間分割数 N, すなわち，空間刻み固定して時間刻み求めた限界の$\Delta t$と$\Delta t$を変えていき，数値計算が可能な最も大きい$\Delta x$の関係を図に示した． を固定して$\Delta x$$\Delta t$$\Delta t$$\Delta x$をを求める．このようにしてを動かした時，この線より大きな$\Delta t$を選ぶと，数値解は発散し計算は不可能となる．

41 バーガース方程式
バーガース方程式の初期値問題を 4 次のルンゲークッタ法，4 次のルンゲークッタ法と積分因子法の組み合わせ法を用いて解いた．このとき図 1 に示すように，衝撃波が完全に立ち上った頃$-0.01$ とした．数値安定性の結果を図 2の時刻 $t=5$ で数値安定性を調べることにした．ただし，「に示す．

図 2 より，4 次のルンゲークッタ法では発散する限界の加は粘性項の 2 階微分の影響から$\Delta x$の約 2 乗に比例し数値安定性はあまり良くなかったが，4 次のルンゲークッタ法と積分因子法の組み合わせ法では$\Delta t$は$\Delta x$に関わらず安定となり，数値安定性は非常に改善された．このため，

231
この場合は他の方法は試みなかった．

42 $KdV$ 方程式

方程式の初期値問題を，4 次のルンゲークッタ法，4 次のルンゲークッタ法と積分因子法の組み合わせ法，線形項と非線形項を段階的に解く方法 (非線形項は 4 次のジェイムソンの方法で解いた．) を用いて解いた．このとき図 3 に示すように，振幅が最も大きく，分散性が効いてくるという $KdV$ 方程式の特徴が出てくる頃の時刻「$-5$ で数値安定性を調べることにした．ここで，$KdV$$\mu=0.01$とした．数値安定性の結果を図 4 に示す．

図 4 より，4 次のルンゲークッタ法では発散する限界の の値は分散項の 3 階微分によりのおよそ 3 乗に比例し，数値安定性は悪かった．積分因子法をルンゲークッタ法との組み合わせで用いるとの限界はのおよそ 2 乗に比例し，かなり改善することが確認できた．さらに，$\Delta t$$\Delta t$$\Delta x$$\Delta x$

線形項と非線形項を段階的に解くと$\Delta t$の限界は$\Delta x$にほぼ比例するまで良くなることが分かった．

ここで，積分因子法とルンゲークッタ法を組み合わせた方法と線形項と非線形項を段階的に解く
方法では図 4 のようにグラフの線が直線ではないので，条件の厳しい値をとるようそれぞれ左から 3 点と 2 点のみを用いて累乗近似曲線の式を求めた．

4.3 KdV-バーガース方程式
粘性項のあるバーガース方程式ではルンゲークッタ法と積分因子法を組み合わせた数値スキームを用いると$\Delta t$

は$\Delta x$に関わらず安定にすることができたが，分散項のある $KdV$ 方程式ではどの数値スキームをしても数値安定となる限界の$\Delta t$

は $O(\Delta t)\sim o(\Delta x)$ よりも安定性が向上することはなかった．では，粘性項と分散項を合わせ持つ KdV- バーガース方程式ではどうなるかを調べた．
KdV- バーガース方程式の初期値問題を，4 次のルンゲ クッタ法，4 次のルンゲークッタ法と-積分因子法の組み合わせ法を用いて解いた．このとき図 5 のように振幅が最も大きく，分散性が効いてくる頃の時刻 $t=5$ で数値安定性を調べることにした．ただし，$v=0.O1,$ $\mu=0.01$ とした．数値安定性の結果を図 6 に示す．

図 6 より，4 次のルンゲークッタ法では $KdV$ 方程式の時と同様に発散する限界の階微分の影響で$\Delta x$$\Delta t$の値は 3の約 3 乗に比例し，数値安定性はあまり良くなかった．4 次のルンゲークッタ法と積分因子法の組み合わせ法ではバーガース方程式の時と同様に$\Delta t$は$\Delta x$に関わらず数値安定になった．このことより，KdV- バーガース方程式は粘性項と分散項を合わせ持っが，粘性項のみを持つバーガース方程式と同様，ルンゲークッタ法と積分因子法を組み合わせた数値スキームを用いると$\Delta t$は$\Delta x$に関わらず数値安定になった．

4.45th $KdV$ 方程式
3 階微分の分散項のある $KdV$ 方程式では線形項と非線形項を段階的に解くと発散しない限界の$\Delta t$の値は$\Delta x$に比例した．5 階微分でもこの方法で数値安定の限界の$\Delta t$が改善するか調べた．

5th $KdV$ 方程式の初期値問題を 1 次，4 次のルンゲークッタ法と積分因子法の組み合わせ法，

232
線形項と非線形項を段階的に解く方法 (非線形項は 4 次のジェイムソンの方法で解いた．) を用いて解いた．このとき図 7 のように振幅が最も大きく，分散性が効いてくるという 5th $KdV$ 方程式の特徴が出てくる頃の時刻 $t=5$ で数値安定性を調べることにした．ただし，$\gamma=0.01$とした．数値安定性の結果を図 8 に示す．

図 8 より，5th $KdV$ 方程式を解くとき，1 次のルンゲークッタ法と積分因子法の組み合わせ法のおよそ 34 乗に比例し，数値安定性はあまり良くなかった．
では発散する限界の の値はここで，このときの結果は図 8 のように線が曲がっているので，条件の厳しい値をとるよう左から 2 点のみを用いて累乗近似曲線の式を求めた．4 次のルンゲークッタ法と積分因子法の組み合はの約 2 乗に比例し，さらに，線形項と非線形項を分けて解く方法ではわせ法では限界の$KdV$ 方程式の時と同様にはにほぼ比例するまで改善した．今回は行っていないが，4 次の$\Delta x$$\Delta t$$\Delta t$$\Delta x$$\Delta t$$\Delta x$

ルンゲークッタ法を用いて解いた場合は，5 階微分の影響により，発散する限界の$\Delta t$は$\Delta x$のおよそ 5 乗に比例すると推測できる．

45 精度について精度を調べるために式 (4), (5) に示す $KdV$ 方程式の保存量である質量 $M$ とエネルギー $E$ の初期$-50$ の値の相対誤差を比較した．時間発展には線形項と非線形項を段階的に解く方法 (非線値と「形項をジェイムソンの方法で解いた．) を用い，アライアジングエラーは 213 則によって除去した．

ジェイムソンの方法の精度を 2 次精度から 10 次精度まで動かすと保存量の相対誤差がどうなるかを調べ，表 1 にまとめた．表 1 より，4 次以上の精度にすると理論的には誤差は減るはずだが，実際には誤差は減らず，精度は改善されなかった．しかし，図 9 に示すように次数をあげると数値安定性はわずかに改善された．ここで，図 9 はある次数において発散する限界の時間刻みが$-C\Delta x$「$\Delta$となる時，横軸に次数，縦軸に表1$C$をとったものである．

各精度における質量とエネルギーの相対誤差

5 まとめ$KdV-$ バーガース方程式はルンゲークッタ法と積分因子法の組み合わせ法
バーガース方程式，$KdV$ 方程式，5th $KdV$ 方程式はこの方法では十分な数値で十分な数値安定性が得られた．一方，安定性は得られず，線形項と非線形項を段階的に解くことで十分な数値安定性が得られた．
線形項と非線形項を段階的に解き，非線形項にジェイムソンの方法を用いることで他の方法を

233用いるときと比べ格段に容易に高次にすることが可能であるが，実際に計算してみると精度はあまり改善しなかった．しかし，高次にすることで数値安定性はわずかに改善した．
6 参考文献[11 J. M. Burgers, “A mathematical model illuminating the theory of turbulence.” Adv. Appl.

Mech. 1, 171-199 (1948).
[2] D. J. Korteweg, G. de Vries, “On the change of form of long waves advancing in arectangular canal and on a new type of long stationary waves.” Philos. Mag. 39, 422-443(1895).
[3] C. Klein, C. Sparber, P. Markowich, “Numerical Study of Oscillatory Regimes in the$Kadomtsev^{-}Petviashvih$ Equation.” Journal of Nonlinear Science 17, 429-470 (2007)[4] R. S. Johnson, “A $non^{-}1inear$ equation incorporating damping and dispersion.” J. Fluid

Mech. 42, 49-60 (1970).
[51 A. Jameson, H. Schmidt, E. Turkel, “Numerical Solutions of the Euler Equations by Finite

Volume Methods Using Runge-Kutta Time Stepping Schemes.” AIAA Paper No. 81-1259(1981).
[6] R. S. Rogallo, ‘An ILLIAC Program for the Numerical Simulation of Homogeneous,

Incompressible Turbulence.” NASA TM $-73203$ (1977).

$0$$\ovalbox{\tt\small REJECT}$

234561

図1バーガース方程式の速度分布 (多-5)

図2$\Delta x-\Delta t$

グラフ (バーガース方程式)

234

1

図3$KdV$ 方程式の速度分布 ( -5)$\ulcorner$$1.E-041.E-031.E-021.E$ $01$ 解く方汐ー$\Delta x$

図4$\Delta x-\Delta t$

グラフ$(KdV$ 方程式$)$1

図 5 $KdV-$ バーガース方程式の速度分布 ( -5)$\ulcorner$

図6$\Delta x-\Delta t$

グラフ ( $KdV-$ バーガース方程式)

235$0$$|$234561

図75th $KdV$ 方程式の速度分布 ( -5)$\ulcorner$

図8$0$

264
次数
図9

次数と数値安定性$\Delta x-\Delta t$

グラフ (5th $KdV$ 方程式)

810