
% Mesinger and Arakawa
% 20110310 (担当者：川畑拓也)



\documentclass[a4j,12pt]{jreport}
\usepackage{Dennou6}
\usepackage{amsmath}
\usepackage{times}
\usepackage{url}
\usepackage{bm}
\pagestyle{Dmyheadings}

\setlength{\parskip}{2ex}

\setcounter{chapter}{10}
\setcounter{section}{1}
\setcounter{subsection}{0}
\setcounter{figure}{0}
\setcounter{table}{0}
\setcounter{page}{1}
\setcounter{equation}{0}

\makeatletter  % 再定義する際に必要. documentclass の後ろに書く必要があ
       % る.
\renewcommand{\chapter}{%
  \thispagestyle{jpl@in}%
  \global\@topnum\z@
  \@afterindenttrue
  \secdef\@chapter\@schapter}
 % chapter の直前で改行しないように設定.
\makeatother

\def\thesection{\arabic{chapter}.\arabic{section}}
\def\theequation{\arabic{section}.\arabic{equation}}
\def\thepage{\arabic{page}} % 数字
%\def\thepage{\roman{page}} % ローマ数字
\def\thetable{\arabic{chapter}.\arabic{section}.\arabic{table}}
\def\thefootnote{\arabic{footnote} )}
\def\thefigure{\arabic{chapter}.\arabic{section}.\arabic{figure}}

\Dtitle{GFDワークノート}{}
\Dauthor{川畑拓也}{}
\Ddate{2011/03/17}
\markright{第11章 アダムズ-バッシュフォース公式}


\begin{document}
\chapter{アダムズ-バッシュフォース公式}
\section{多段階公式}

オイラー法, ルンゲ-クッタ型公式は $t = (n+1) \Delta t$ の従属変数の値を
予報するのに, $t = n\Delta t$ の従属変数の値を用いてきた (1 段階公式). 
これに対し, アダムズ-バッシュフォース公式は $t = (n+1) \Delta t$ の従属変数
の値を予報するのに, $t = n \Delta t$, $(n-1) \Delta t$, $(n-2) \Delta t$, 
$\cdots$ の従属変数の値を用いる (多段階公式). 

解く式は, 
\begin{equation}
  \DD{x}{t} = f(x,t) \Deqlab{eq1}
\end{equation}
とする. 時間刻みを $h$ として,
\begin{equation}
  x(t_n + h) 
\end{equation}
を, $t = t_n$ のまわりで展開する. 展開式は,
\begin{align}
  x(t_n + h) &= x(t_n) + \left( \DD{x}{t} \right)_n h
                       + \frac{1}{2} \left( \DD[2]{x}{t} \right)_n h^2 + O(h^3) \notag \\
             &= x(t_n) + h f + \frac{1}{2} h^2 \left( \DD{f}{t} \right)_n + O(h^3)
  \Deqlab{tenkai}
\end{align}
ルンゲ-クッタ型公式は, $f_n$ をいろいろな $t_n$ での $(x_n, t_n)$ で評価してその重み
つき平均で, \Deqref{tenkai}を近似する. アダムズ-バッシュフォース公式では, $f^{\prime} = \DD{f}{t}$, 
$f^{\prime \prime}$ を $t_n$ よりも前の時刻での $f_n$ の値を用いて, 
過去に計算された $f$ の値を用いて計算し, 展開式 \Deqref{tenkai}を近似する. 

\section{2 位のアダムズ-バッシュフォース公式}
\begin{equation}
  f^{\prime} = \DD{f}{t} \approx \frac{f_n - f_{n-1}}{n \Delta t}
\end{equation}
と近似する. $f_{n-1}$ は $x_{n-1}$ が $O(h^2)$ で求められているとすると,
\begin{align}
  f_{n-1} &= f(x_{n-1}, t_{n-1}) \notag \\
          &= f \left( x(t_n-h) + O(h^2), t_n -h \right) \notag \\
          &= f \left( x(t_n-h), t_n-h \right) + O(h^2). \Deqlab{2.5} 
\end{align}
\Deqref{2.5} を $t_n$ のまわりで展開すると, 
\begin{equation*}
  f_{n-1} = f(x_n, t_n) - h f^{\prime}(x_n, t_n) + O(h^2).
\end{equation*}
ゆえに,
\begin{equation}
  f^{\prime} = \frac{f_n - f_{n-1}}{h} + O(h). \Deqlab{2.6}
\end{equation}
\Deqref{2.6} を展開式 \Deqref{tenkai} に代入すると,
\begin{align}
  x(t_n + h)
  &= x(t_n) + h f_n + \frac{1}{2} h^2 \left[ \frac{f_n - f_{n-1}}{h} + O(h) \right] + O(h^3) \notag \\
  &= x(t_n) + h f_n + \frac{1}{2} h \left[ f_n - f_{n-1} + O(h^2) \right] + O(h^3) \notag \\
  &= x(t_n) + h \left( \frac{3}{2} f_n - \frac{1}{2} f_{n-1} \right) + O(h^3). \Deqlab{2.7}
\end{align}
ここで,
\begin{align*}
  f_n &= f(x_n, t_n) \\
  f_{n-1} &= f(x_{n-1}, t_{n-1})
\end{align*}
である. 

\section{3 位のアダムズ-バッシュフォース公式}
2 位の時と同様にして, $x(t_n + h)$ を 3 次の項まで展開する.
\begin{equation}
  x(t_n + h) = x(t_n) + hf_n + \frac{1}{2} h^2 f^{\prime} + \frac{1}{6}h^3 f^{\prime \prime} + O(h^4). \Deqlab{tenkai2} 
\end{equation}
同様に,
\begin{align*}
  f(x_{n-1}, t_{n-1}) &= f(x(t_n - h), t_n -h ) + O(h^3) \\
  &= f_{n-1} + O(h^3).
\end{align*}
ゆえに,
\begin{align}
  f_{n-1} &= f_n -h f_n^{\prime} + \frac{1}{2} h^2 f_n^{\prime \prime} + O(h^3) \Deqlab{1} \\
  f_{n-2} &= f_n - 2h f_n^{\prime} + \frac{1}{2} (2h)^2 f_n^{\prime \prime} + O(h^3). \Deqlab{2}
\end{align}
\Deqref{1} と \Deqref{2} より,
\begin{equation*}
  4f_{n-1} - f_{n-2} = 3f_n - 2h f_n^{\prime} + O(h^3).
\end{equation*}
従って, 
\begin{equation}
  f_n^{\prime} = \frac{3f_n -4f_{n-1} + f_{n-2}}{2h} + O(h^2). \Deqlab{3}
\end{equation}
同様に \Deqref{1}, \Deqref{2} より,
\begin{equation*}
  2f_{n-1} - f_{n-2} = f_n - h^2 f_n^{\prime \prime} + O(h^3)
\end{equation*}
であるから,
\begin{equation}
  f_n^{\prime \prime} = \frac{f_n - 2f_{n-1} + f_{n-2} }{h^2} + O(h^3). \Deqlab{4}
\end{equation}
\Deqref{3}, \Deqref{4} を \Deqref{tenkai2} へ代入すると,
\begin{align}
  x_{n+1}
  &= x_n + hf_n + \frac{1}{2} h^2 \left( \frac{3f_n - 4 f_{n-1} + f_{n-2}}{2h} \right)
  + \frac{1}{6} h^3 \left( \frac{f_n -2 f_{n-2} + f_{n-2}}{h^2} \right) \notag\\
  &= x_n + h \left( \frac{23}{12}f_n - \frac{16}{12} f_{n-1} + \frac{5}{12} \right) \notag \\
  &= x_n + \frac{h}{12} ( 23 f_n - 16 f_{n-1} + 5 f_{n-2})
\end{align}
となる. これが 3 位のアダムズ-バッシュフォースの公式である.


\section{アダムズ-バッシュフォース公式の特徴}
\begin{description}
  \item[長所]\mbox{}
    \begin{itemize}
      \item 何位の公式でも (ルンゲ-クッタ型公式に比べ) 簡単に作れる. なぜならば
            ステップ数とオーダーが必ず一致するからである.
      \item 1 ステップあたりの関数評価が 1 回で済む
    \end{itemize}
  \item[短所]\mbox{}
    \begin{itemize}
      \item 高次になるにつれて安定領域が狭くなる.
      \item 計算開始時には過去の関数値 $f_k$ を使えないので, 別の方法であらかじめ $f_k$ 
            を求めておく必要がある.
    \end{itemize}
\end{description}

\section{予測子・修正子法}

アダムズ-バッシュフォース公式では, $f^{\prime}$, $f^{f^{\prime \prime}}$ を求める
のに, 過去の時刻の $f$ の値を用いた. これに対し, 予測子・修正子法では
 $f^{\prime}$, $f^{\prime \prime}$ を求めるのに, 予報時刻の $f$ の値を用いる. 
例えば 2 位の公式では,
\begin{equation*}
  f(x_{n+1}, t_{n+1}) = f(x_n, t_n) + h f^{\prime} (x_n, t_n) + O(h^2)
\end{equation*}
より,
\begin{equation}
  f_n^{\prime} = \frac{f_{n+1} -f_n}{h}. \Deqlab{fprime}
\end{equation}
\Deqref{fprime} を \Deqref{tenkai} に代入すると,
\begin{align}
  x(t_n + h) 
  &= x(t_n) + hf_n + \frac{1}{2}h^2 \frac{f_{n+1} - f_n}{h} \notag \\
  &= x(t_n) + hf_n + \frac{1}{2}hf_{n+1} - \frac{1}{2}hf_n \notag \\
  &= x(t_n) + \frac{1}{2}h(f_{n+1} + f_n).
\end{align}
ゆえに,
\begin{equation*}
  x_{n+1} = x_n + \frac{h}{2} (f_{n+1}+ f_n).
\end{equation*}
ただし,
\begin{equation}
  \begin{cases}
    f_{n+1} &= f(x^{\ast}, t^{\ast}) \\
    x^{\ast} &= x_n + h \left( \frac{3}{2} f_n - \frac{1}{2} f_{n-1} \right) \\
    t^{\ast} &= t_n + h
  \end{cases}
\end{equation}
となる.

\section{アダムズ-バッシュフォース公式の安定性}
振動方程式
\begin{equation}
  \DD{x}{t} = \omega x \ \qquad \ (f \equiv \omega x)
\end{equation}
にアダムズ-バッシュフォース公式を用いた場合の安定性を調べる.
\subsection{2 位の公式の場合}
\begin{align*}
  x_{n+1} &= x_n + h \left( \frac{3}{2} f_n - \frac{1}{2} f_{n-1} \right) \\
  &= x_n + h \left( \frac{3}{2} \omega x_n - \frac{1}{2} \omega x_{n-1} \right).
\end{align*}
増幅係数 $\lambda$ を,
\begin{equation}
  \lambda \equiv \frac{x_{n+1}}{x_n} \Deqlab{lambda}
\end{equation}
と定義する. このとき,
\begin{equation*}
  \lambda = 1 + h \left( \frac{3}{2} \omega - \frac{1}{2} \frac{\omega}{\lambda} \right)
\end{equation*}
であるから,
\begin{equation*}
  \lambda^2 - \left( 1 + \frac{3}{2} h \omega \right) + \frac{1}{2} h \omega =0.
\end{equation*}
$h \omega \equiv p$ とおき, 上式を $\lambda$ について解くと,
\begin{equation*}
  \lambda^2 - \left( 1 + \frac{3}{2} p \right) \lambda + \frac{1}{2} p =0
\end{equation*}
より,
\begin{equation}
  \lambda = \frac{1}{2} 
  \left( 1 + \frac{3}{2} p \pm \Dsqrt{1 + p + \frac{9}{4}p^2} \right). \Deqlab{lambda_kai}
\end{equation}
$p \ll 1$ とすると, \Deqref{lambda_kai} は,
\begin{align*}
  \Dsqrt{1+p+\frac{9}{4}p^2}
  &= 1 + \frac{1}{2} \left( p+ \frac{9}{4}p^2 \right) - \frac{1}{8} \left( p+ \frac{9}{4}p^2 \right)^2 + \cdots \\
  &= 1+ \frac{1}{2} p + p^2 + O(p^3).
\end{align*}
よって, $\lambda$ の 2 つの解の $p$ の 2 次までの近似式は,
\begin{align}
  \lambda_+ &= 1+p+ \frac{1}{2}p^2 +O(p^3), \\
  \lambda_- &= \frac{1}{2}p - \frac{1}{2}p^2 + O(h^3).
\end{align}
$x_{n+1} = \lambda x_n$ より, $p \to 0$ の極限で,
\begin{align*}
  \lambda_+ x_n &\to x_{n+1}, \\
  \lambda_- x_n &\to 0.
\end{align*}
これより, $\lambda_+$ に対応する数値解は物理モード, $\lambda_-$ に対応する数値解は
計算モードである. 安定性の条件は $\lvert \lambda \rvert \le 1$, $p \ll 1$ として, 
$p$ の 1 次の項まで考えると,
\begin{align*}
  \lvert \lambda_+ \rvert &= \Dsqrt{(1+p_{\mathrm{re}})^2 + p_{\mathrm{im}}^2} > 1, \\
  \lvert \lambda_- \rvert &= \frac{1}{2} \lvert p \rvert < 1.
\end{align*}
ここで, $p_{\mathrm{re}} = \mathrm{Re[p]}$, $p_{\mathrm{im}} = \mathrm{Im[p]}$ 
である. 
よって, 2 位公式の場合, 物理モードは不安定である. 一方計算モードは安定となる.
安定領域の境界は $\lvert \lambda \rvert =1$ なので, 
\begin{equation*}
  p = \frac{\lambda^2 - \lambda}{\frac{3}{2}\lambda - \frac{1}{2}}
\end{equation*}
に $\lambda = e^{i \theta}$ を代入して, $0 \le \theta \le 2 \pi$ とすれば得られる.

\subsection{2 位の予測子・修正子法の安定性}
\begin{align*}
  x_{n+1} &= x_n + \frac{h}{2} (\omega x^{\ast} + \omega x_n), \\
  x^{\ast} &= x_n + h \left( \frac{3}{2} \omega x_n - \frac{1}{2}\omega x_{n-1} \right)
\end{align*}
より,
\begin{equation}
  x_{n+1} = \left( 1 + h\omega + \frac{4}{3} (h\omega)^2 \right) x_n - \frac{1}{4} (h \omega)^2 x_{n-1}.
\end{equation}
再び増幅係数 $\lambda$ を,
\begin{equation*}
  \lambda \equiv \frac{x_{n+1}}{x_n}
\end{equation*}
とおくと,
\begin{equation*}
  \lambda^2 - \left( 1 + p + \frac{3}{4} p^2 \right) \lambda + \frac{1}{4} p^2 = 0.
\end{equation*}
ただし, $p = h \omega$ である. これを $\lambda$ について解くと,
\begin{equation*}
  \lambda = \frac{1}{2} 
  \left(
  1 + p + \frac{3}{4} p^2 \pm \Dsqrt{1+2p+\frac{3}{2}p^2+\frac{3}{2}p^3 + \frac{9}{16}p^4}
  \right).
\end{equation*}
$p \ll 1$ として,
\begin{align*}
  \Dsqrt{1+2p+\frac{3}{2}p^2+\frac{3}{2}p^3 + \frac{9}{16}p^4}
  &= 1 + \frac{1}{2} \left( 2p+\frac{3}{2}p^2 + \frac{3}{2}p^3 + \frac{9}{16}p^4 \right) \\
  &-\frac{1}{4} \left( 2p+\frac{3}{2}p^2 + \frac{3}{2}p^3 + \frac{9}{16}p^4 \right)^2 + O(p^9) \\
  &= 1 + p - \frac{1}{4} p^2 + O(p^3).
\end{align*}
ゆえに,
\begin{align}
  \lambda_+ &= 1+ p -\frac{1}{4}p^2 +O(p^3), \\
  \lambda_- &= \frac{1}{2} p^2 + O(p^3).
\end{align}
2 位のアダムズ-バッシュフォース公式と比較すると, 物理モードの振る舞いはほぼ同じであり, 
計算モードの振る舞いは, 予測子・修正子法の方が早く減衰する ($\lvert p \rvert < 1$ のとき).
\end{document}
