% % 表題 SPMODEL サンプルプログラムドキュメント % 1 次元 KdV 方程式モデル % % 履歴 2004/02/25 小高 正嗣 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \documentclass[a4j,12pt]{jarticle} % % プリアンブル % \usepackage{Dennou6} % dennou-sty-6 を使用 \Dtitle[1 次元 KdV 方程式モデル] % ヘッダ部タイトル, モデル名と同じ { {\large SPMODEL サンプルプログラム} \\ % ここは変更しない 1 次元 KdV 方程式モデル \\ % モデル名 {\large kdv1.f90, kdv2.f90} % プログラム名 } % \Dauthor[小高 正嗣]{小高 正嗣} % 著者 \Ddate{2004 年 2 月 25 日} % 日付 \Dnoparindent % 段落の字下げをしない % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % 本文開始 % % 基本的な章節構成は以下の通り. % % 1. はじめに % 2. 支配方程式系 % 3. 離散化 % 4. 利用モジュールとその他の設定 % 5. 数値実験 % 6. 参考文献 % -. 謝辞 % % 必要に応じて構成を変更すること. % \begin{document} \maketitle % タイトルの作成 \tableofcontents % 目次の作成 \pagebreak %---------------------------------------------------------------------- \Dparskip %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{概要} SPMODEL サンプルプログラム『kfv1.f90』および『kfv2.f90』に用いられてい る基礎方程式と境界条件, および, このプログラムを用いた数値実験の方法に ついて解説する. 基礎方程式は, 1 次元の KdV 方程式である.計算はスペクト ル法を用いて行い, 展開はフーリエ級数を用いて行う. スペクトル変換 と逆変換および微分演算には, SPMODEL ライブラリ(spml)を用いている. 時間 積分法は『kfv1.f90』では leapfrog スキーム, 『kfv2.f90』では Adams-Bashforth スキームを用いる. 数値実験では振幅の異なる 2 つの孤立 波解を初期値に与え, それらの時間発展を計算する. {\bf プログラム名} \\ % プログラム名 {\footnotesize kdv1.f90, kdv2.f90 } {\bf プログラム取得元}\\ % プログラム取得先 {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash 1d-cyclic-e\slash kdv\slash sample\slash SIGEN.htm } {\bf SPMODEL サンプルプログラム目次}\\ % サンプルプログラム目次(変更しない) {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash sample.htm } {\bf SPMODEL の使い方}\\ % SPMODEL の使い方(変更しない) {\footnotesize http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\/ } %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{支配方程式系} ここでは支配方程式系と境界条件を記す. \subsection{支配方程式系} 支配方程式系は 1 次元の KdV 方程式である: \begin{eqnarray} \DP{\zeta}{t} + \zeta\DP{\zeta}{x} + \DP[3]{\zeta}{x} = 0. \end{eqnarray} 各記号は以下の量をあらわす (以下のように表にまとめる). \begin{table}[h] \begin{center} \begin{tabular}{c l l } \hline 記号 &\qquad\qquad& 変数/物理定数 \\ \hline $x$ && 水平位置 \\ $t$ && 時間 \\ $\zeta$ && 振幅 \\ \hline \end{tabular} \end{center} \caption{変数, 物理定数の定義} % \Dtablab{} \end{table} \subsection{境界条件} 水平方向の境界条件は周期境界条件である. すなわち, 水平計算領域を $x_{m}$ とすると, \[ \zeta (x+x_{m}) = \zeta (x), \] である. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{離散化} この節では方程式系の空間離散化および使用した時間積分法について説明し, プログラム内で実際に用いられている方程式を記述する. \subsection{空間離散化} 支配方程式の離散表現は以下のようになる. \begin{eqnarray} \DP{\zeta_{i}(t)}{t} + \zeta _{i}(t)\DP{\zeta_{i}(t)}{x} + \DP[3]{\zeta_{i}(t)}{x} = 0. \Deqlab{空間離散化した式} \end{eqnarray} ここで下付き添字 $i$ は水平格子点 $x_{i}$ 上での値を示す. \subsection{水平方向のスペクトル表現} 空間離散化した支配方程式\Deqref{空間離散化した式}をスペクトル法を用い て表現する. 空間方向のスペクトル展開はフーリエ級数を用いて行う. ただし 非線形項を扱う場合は, 先に格子点上での非線形項の値を計算し, その値のス ペクトルを求める方法(変換法)を用いる. 以下では $k$ を $x$ 方向波数, $K$ を切断波数, $I$ を $x$ 方向格子点数とする. $\zeta_{i}(t)$ はスペクトル逆変換によって以下のように展開される: \begin{equation} \zeta_{i}(t) = \sum _{k=0}^{K} \exp\left(\frac{2\pi i k x_{i}}{x_{m}}\right)\hat \zeta_{k}. \Deqlab{スペクトル逆変換} \end{equation} スペクトル係数 $\hat \zeta_{k}$ は以下のスペクトル変換によって与えられる: \begin{equation} \hat \zeta_{k} = \frac{1}{I} \sum _{i=0}^{I-1}\exp\left(-\frac{2\pi i k x_{i}}{x_{m}}\right)\zeta_{i}. \Deqlab{スペクトル変換} \end{equation} \Deqref{スペクトル逆変換}, \Deqref{スペクトル変換}より, \Deqref{空間離散化した式}のスペクトル表現は以下のようになる. \begin{equation} \DP{\hat \zeta_{k}}{t} = - \frac{1}{I} \sum _{i=0}^{I-1}\exp\left(-\frac{2\pi i x_{i}}{x_{m}}\right) \zeta _{i}\left(\DP{\zeta_{i}}{x}\right)_{i} \\ - \left(\frac{2\pi i k}{x_{m}}\right)^{3}\hat \zeta_{k}. \Deqlab{支配方程式のスペクトル表現} \end{equation} \subsection{時間積分} ここでは時間積分法について記述し, プログラム内で実際に用いられている方 程式を記述する. 以下では $\Delta t$ を時間格子間隔, 時刻 $\tau \Delta t$ における $\hat \zeta_{k}$ の値を $\hat \zeta_{k}^{\tau}$ 等と表す. \subsubsection{leapfrog スキーム} leapfrog スキームを用いて\Deqref{支配方程式のスペクトル表現}式を時間方 向に離散化すると, 以下のように表される. \begin{equation} \hat \zeta_{k}^{\tau +1} = \hat \zeta_{k}^{\tau -1} + 2\Delta t \hat F_{k}^{\tau}. \end{equation} ここで, \begin{equation} \hat F_{k}^{\tau} = - \frac{1}{I} \sum _{i=0}^{I-1}\exp\left(-\frac{2\pi i x_{i}}{x_{m}}\right) \zeta _{i}^{\tau}\left(\DP{\zeta_{i}}{x}\right)_{i}^{\tau} - \left(\frac{2\pi i k}{x_{m}}\right)^{3}\hat \zeta_{k}^{\tau} \Deqlab{時間積分右辺} \end{equation} である. 1 ステップ目を計算する場合は Euler スキームを用いる. \subsubsection{Adams-Bashforth スキーム} Adams-Bashforth スキームを用いて\Deqref{支配方程式のスペクトル表現}式 を時間方向に離散化すると, 以下のように表される. \begin{equation} \hat \zeta_{k}^{\tau +1} = \hat \zeta_{k}^{\tau } + \frac{\Delta t}{2} \left( 3\hat F_{k}^{\tau} - \hat F_{k}^{\tau -1}\right). \end{equation} $F_{k}^{\tau}$ の表現は\Deqref{時間積分右辺}式と同じである. 1 ステップ目を計算する場合は Euler スキームを用いる. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{使用モジュールとその他の設定} スペクトル変換と逆変換, 微分演算は SPMODEL ライブラリ (spml) の ae\_module に含まれる関数を用いて行う. spml が下位で使用する ISPACK の 仕様から, 格子点数 $I$ は偶数で, かつ $I/2=2^{a}3^{b}5^{c}$ ($a, b, c$ は 0 または整数) でなければならない. 非線形項の計算によって生じるエリ アジングを防ぐため, 格子点数 $I$ と切断波数 $K$ は $I>3K$ を満たすよう に与える. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section{数値実験} 初期に振幅の異なる 2 つの孤立波を与え, それらの時間発展を計算する. 初 期条件は以下のように与える(\Dfigref{init}). \begin{equation} \zeta = U_{1}\mbox{sech}^{2}\left(\frac{x-x_{1}}{\sqrt{12/U_{1}}}\right) + U_{2}\mbox{sech}^{2}\left(\frac{x-x_{2}}{\sqrt{12/U_{2}}}\right). \Deqlab{初期条件} \end{equation} $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 スキームを用いる場合)では $\Delta t=10^{-6}$ sec, 『kdv2.f90』 (Adams-Bashforth スキームを用いる場合)で は$\Delta t=2\times 10^{-7}$ sec である. 計算ステップ数はそれぞれ 5,000と 250,000 である. \Dfigref{kdv1},\Dfigref{kdv2}に計算結果を示す. 振幅が大きいほど位相速 度も大きいこと, 振幅の大きい波が振幅の小さい波を追い越して行く様子がわ かる. \begin{figure}[b] \begin{center} \Depsf[][90mm]{./figs/init.ps} \caption{初期値の分布.} \Dfiglab{init} \end{center} \end{figure} \begin{figure}[p] \begin{center} \Depsf[][100mm]{./figs/kdv1.ps} \caption{Leapfrog スキームを用いた場合の数値計算の結果. $x$-$t$ 断面図.} \Dfiglab{kdv1} \Depsf[][100mm]{./figs/kdv2.ps} \caption{Adams-Bashforth + Crank Nicholson スキームを用いた場合の数値計算の結果. $x$-$t$ 断面図.} \Dfiglab{kdv2} \end{center} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \newpage \section{参考文献} \begin{description} \item Witham, G. B., 1974: {\it Linear and nonlinear waves}, John Wiley \& Sons, 636pp. \item 巽 友正, 1982: 流体力学, 培風館. \item 竹広真一, 石岡圭一, 豊田英司, 石渡正樹, 林祥介, SPMODEL 開発グループ, 2002: 階層的地球流体力学スペクトルモデル集 (SPMODEL),\\ http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash , 地球流体電脳倶楽部. \end{description} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newpage % 改ページ \section*{謝辞} % 謝辞. このまま利用する. \markright{謝辞} % ヘッダに表示 \addcontentsline{toc}{section}{謝辞} % 目次に表示 本資源は, 地球流体電脳倶楽部のインターネット上での学術知識の集積と活用 の実験の一環として \begin{center} http:\slash \slash www.gfd-dennou.org\slash arch\slash spmodel\slash \end{center} において公開されているものである (\copyright 地球流体電脳倶楽部スペクトルモデルプロジェクト spmodel@gfd-dennou.org 2002. ). 本資源は, 著作者の諸権利に抵触しない(迷惑をかけない)限りにおいて自由に 利用していただいて構わない. なお, 利用する際には今一度自ら内容を確かめ ることをお願いする(無保証無責任原則). 本資源に含まれる元資源提供者(図等の版元等を含む)からは, 直接的な形で の WEB 上での著作権または使用許諾を得ていない場合があるが, 勝手ながら, 「未来の教育」のための実験という学術目的であることをご理解いただけるも のと信じ, 学術標準の引用手順を守ることで諸手続きを略させていただいて いる. 本資源の利用者には, この点を理解の上, 注意して扱っていただける ようお願いする. 万一, 不都合のある場合には \begin{center} spmodel@gfd-dennou.org \end{center} まで連絡していただければ幸いである. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \end{document}