\underline{1 次元 KdV 方程式モデル目次 1}

SPMODEL サンプルプログラム

1 次元 KdV 方程式モデル

kdv1.f90, kdv2.f90

小高正嗣

2004 年 2 月 25 日

目次

1 概要 2 2 支配方程式系 3

2.1 支配方程式系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3 離散化 4 3.1 空間離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.2 水平方向のスペクトル表現 . . . . . . . . . . . . . . . . . . . . . . . 4 3.3 時間積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.3.1 leapfrog スキーム . . . . . . . . . . . . . . . . . . . . . . . . 5 3.3.2 Adams-Bashforth スキーム . . . . . . . . . . . . . . . . . . 5

4 使用モジュールとその他の設定 6 5 数値実験 7 6 参考文献 9 謝辞 10

\underline{1 次元 KdV 方程式モデル  1 概要  2}

1 概要

SPMODEL サンプルプログラム『kfv1.f90』および『kfv2.f90』に用いられている基礎方程式と境界条件, および, このプログラムを用いた数値実験の方法について解説する. 基礎方程式は, 1 次元の KdV 方程式である. 計算はスペクトル法を用いて行い, 展開はフーリエ級数を用いて行う. スペクトル変換と逆変換および微分演算には, SPMODEL ライブラリ (spml) を用いている. 時間積分法は『kfv1.f90』では leapfrog スキーム, 『kfv2.f90』では Adams-Bashforth スキームを用いる. 数値実験では振幅の異なる 2 つの孤立波解を初期値に与え, それらの時間発展を計算する.

プログラム名

kdv1.f90, kdv2.f90

プログラム取得元

http://www.gfd-dennou.org/arch/spmodel/1d-cyclic-e/kdv/sample/SIGEN.htm

SPMODEL サンプルプログラム目次

http://www.gfd-dennou.org/arch/spmodel/sample.htm

SPMODEL の使い方

http://www.gfd-dennou.org/arch/spmodel

\underline{1 次元 KdV 方程式モデル  2 支配方程式系  3}

2 支配方程式系

ここでは支配方程式系と境界条件を記す.

2.1 支配方程式系

支配方程式系は 1 次元の KdV 方程式である:

\frac{∂ζ}{∂t} + ζ \frac{∂ζ}{∂x} + \frac{∂^{3}ζ}{∂x^{3}} = 0.  (1)

各記号は以下の量をあらわす (以下のように表にまとめる).

\begin{matrix}記号 & 変数/物理定数 \\ x & 水平位置 \\ t & 時間 \\ ζ & 振幅 \end{matrix}表 1: 変数, 物理定数の定義

2.2 境界条件

水平方向の境界条件は周期境界条件である. すなわち, 水平計算領域を x_{m} とすると,

ζ(x + x_{m}) = ζ(x),

である.

\underline{1 次元 KdV 方程式モデル  3 離散化  4}

3 離散化

この節では方程式系の空間離散化および使用した時間積分法について説明し, プログラム内で実際に用いられている方程式を記述する.

3.1 空間離散化

支配方程式の離散表現は以下のようになる.

\frac{∂ζ_{i}(t)}{∂t} + ζ_{i}(t)\frac{∂ζ_{i}(t)}{∂x} + \frac{∂^{3}ζ_{i}(t)}{∂x^{3}} = 0.  (2)

ここで下付き添字 i は水平格子点 x_{i} 上での値を示す.

3.2 水平方向のスペクトル表現

空間離散化した支配方程式 (2) をスペクトル法を用いて表現する. 空間方向のスペクトル展開はフーリエ級数を用いて行う. ただし非線形項を扱う場合は, 先に格子点上での非線形項の値を計算し, その値のスペクトルを求める方法 (変換法) を用いる. 以下では k を x 方向波数, K を切断波数, I を x 方向格子点数とする.

ζ_{i}(t) はスペクトル逆変換によって以下のように展開される:

ζ_{i}(t) = ∑_{k=0}^{K} exp (\frac{2πikx_{i}}{x_{m}} ) ζ̂_{k}.  (3)

スペクトル係数 ζ̂_{k} は以下のスペクトル変換によって与えられる:

ζ̂_{k} = \frac{1}{I} ∑_{i=0}^{I−1} exp (−\frac{2πikx_{i}}{x_{m}} ) ζ_{i}.  (4)

(3), (4) より, (2) のスペクトル表現は以下のようになる.

\frac{∂ζ̂_{k}}{∂t} = −\frac{1}{I} ∑_{i=0}^{I−1} exp (−\frac{2πix_{i}}{x_{m}} ) ζ_{i} (\frac{∂ζ_{i}}{∂x} )_{i} − (\frac{2πik}{x_{m}} )^{3} ζ̂_{k}.  (5)

3.3 時間積分

ここでは時間積分法について記述し, プログラム内で実際に用いられている方程式を記述する. 以下では Δt を時間格子間隔, 時刻 τΔt における ζ̂_{k} の値を ζ̂_{k}^{τ} 等と表す.

3.3.1 leapfrog スキーム

leapfrog スキームを用いて (5) 式を時間方向に離散化すると, 以下のように表される.

ζ̂_{k}^{τ+1} = ζ̂_{k}^{τ−1} + 2ΔtF̂_{k}^{τ}.  (6)

ここで,

F̂_{k}^{τ} = −\frac{1}{I} ∑_{i=0}^{I−1} exp (−\frac{2πix_{i}}{x_{m}} ) ζ_{i}^{τ} (\frac{∂ζ_{i}}{∂x} )_{i}^{τ} − (\frac{2πik}{x_{m}} )^{3} ζ̂_{k}^{τ}  (7)

である.

1 ステップ目を計算する場合は Euler スキームを用いる.

3.3.2 Adams-Bashforth スキーム

Adams-Bashforth スキームを用いて (5) 式を時間方向に離散化すると, 以下のように表される.

ζ̂_{k}^{τ+1} = ζ̂_{k}^{τ} + \frac{Δt}{2} (3F̂_{k}^{τ} − F̂_{k}^{τ−1}) .  (8)

F_{k}^{τ} の表現は (7) 式と同じである.

1 ステップ目を計算する場合は Euler スキームを用いる.

\underline{1 次元 KdV 方程式モデル  4 使用モジュールとその他の設定  6}

4 使用モジュールとその他の設定

スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の ae module に含まれる関数を用いて行う. spml が下位で使用する ISPACK の仕様から, 格子点数 I は偶数で, かつ I/2 = 2^{a}3^{b}5^{c} (a, b, c は 0 または整数) でなければならない. 非線形項の計算によって生じるエリアジングを防ぐため, 格子点数 I と切断波数 K は I > 3K を満たすように与える.

\underline{1 次元 KdV 方程式モデル  5 数値実験  7}

5 数値実験

初期に振幅の異なる 2 つの孤立波を与え, それらの時間発展を計算する. 初期条件は以下のように与える (図 1).

ζ = U_{1}sech^{2} (\frac{x − x_{1}}{\sqrt{12/U_{1}}} ) + U_{2}sech^{2} (\frac{x − x_{2}}{\sqrt{12/U_{2}}} ) .  (9)

U_{1} = 720, U_{2} = 1440, x_{1} = x_{m}/2, x_{2} = x_{m}/4 とする. 水平計算領域の大きさは x_{m} = 3 とする.

格子点数 I と切断波数 K はそれぞれ I = 128, K = 42 とする. 時間格子間隔は, 『kdv1.f90』 (leapfrog スキームを用いる場合) では Δt = 10^{−6} sec, 『kdv2.f90』 (Adams-Bashforth スキームを用いる場合) では Δt = 2 × 10^{−7} sec である. 計算ステップ数はそれぞれ 5,000 と 250,000 である.

図 2, 図 3 に計算結果を示す. 振幅が大きいほど位相速度も大きいこと, 振幅の大きい波が振幅の小さい波を追い越して行く様子がわかる.

図 1: 初期値の分布.

\underline{1 次元 KdV 方程式モデル  5 数値実験  8}

図 2: Leapfrog スキームを用いた場合の数値計算の結果. x-t 断面図. 図 3: Adams-Bashforth + Crank Nicholson スキームを用いた場合の数値計算の結果. x-t 断面図.

\overline{kdv.tex 2004 年 2 月 25 日 (小高正嗣)}

\underline{1 次元 KdV 方程式モデル  6 参考文献  9}

6 参考文献

Witham, G. B., 1974: Linear and nonlinear waves, John Wiley & Sons, 636pp. 巽友正, 1982: 流体力学, 培風館. 竹広真一, 石岡圭一, 豊田英司, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2002:

階層的地球流体力学スペクトルモデル集 (SPMODEL), http://www.gfd-dennou.org/arch/spmodel/, 地球流体電脳倶楽部.

\underline{1 次元 KdV 方程式モデル謝辞 10}

謝辞

本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用の実験の一環として

http://www.gfd-dennou.org/arch/spmodel/

において公開されているものである (\begin{matrix}c \\ ◯ \end{matrix}地球流体電脳倶楽部スペクトルモデルプロジェクト spmodel@gfd-dennou.org 2002. ). 本資源は, 著作者の諸権利に抵触しない (迷惑をかけない) 限りにおいて自由に利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめることをお願いする (無保証無責任原則).

本資源に含まれる元資源提供者(図等の版元等を含む)からは, 直接的な形での WEB 上での著作権または使用許諾を得ていない場合があるが, 勝手ながら, 「未来の教育」のための実験という学術目的であることをご理解いただけるものと信じ, 学術標準の引用手順を守ることで諸手続きを略させていただいている. 本資源の利用者には, この点を理解の上, 注意して扱っていただけるようお願いする. 万一, 不都合のある場合には

spmodel@gfd-dennou.org

まで連絡していただければ幸いである.