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



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

\setlength{\parskip}{2ex}

\setcounter{chapter}{9}
\setcounter{section}{1}
\setcounter{subsection}{1}
\setcounter{figure}{0}
\setcounter{table}{0}
\setcounter{page}{5}
\setcounter{equation}{7}

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

\def\thesecion{\arabic{chapter}.\arabic{section}}
\def\thesection{\arabic{section}}
\def\theequation{\arabic{chapter}.\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{2010/10/21}
\markright{第9章 偏微分方程式の数値解法の基礎 2}


\begin{document}
\subsection{3 段階 (3-level) スキーム}
$U^{n+1}$ を求める式に 3 つの$U^i$が現れるスキーム. ただし, 1 ステップ目の計算 ($U^0$ から $U^1$ を求めるとき) には使えない. 一般的な表現は,
\begin{equation}
U^{n+1} = U^{n-1} + \int_{(n-1)\Delta t}^{(n+1)\Delta t}f(U,t) dt
\end{equation}

{\bf B-1. リープフロッグスキーム (Leapfrog scheme)}

\begin{equation}
U^{n+1} = U^{n-1} + 2 \Delta t f^n
\end{equation}
打ち切り誤差は,
\begin{align*}
U^{n+1} &= U^n + \left. \DD{U}{t} \right|_n \Delta t + \frac{1}{2} \left. \DD[2]{U}{t} \right|_n (\Delta t)^2 + \frac{1}{3!} \left. \DD[3]{U}{t} \right|_n (\Delta t)^3 + O(\Delta t^4),\\
U^{n-1} &= U^n - \left. \DD{U}{t} \right|_n \Delta t + \frac{1}{2} \left. \DD[2]{U}{t} \right|_n (\Delta t)^2 - \frac{1}{3!} \left. \DD[3]{U}{t} \right|_n (\Delta t)^3 + O(\Delta t^4)
\end{align*}
より,
\begin{align*}
\varepsilon_{leap} &= \frac{U^{n+1}-U^{n-1}}{2\Delta t} - f^n\\
                   &= \left. \DD{U}{t} \right|_n + \frac{1}{3} \left. \DD[3]{U}{t} \right|_n + (\Delta t)^2 + O(\Delta t^5) -f^n\\
                   &= O(\Delta t^2).
\end{align*}

{\bf B-2. アダムス-バッシュフォーススキーム (Adams-Bashforth scheme)}

オリジナルのアダムス-バッシュフォーススキームは 4 次精度であるが, ここでは 2 次精度のスキームを紹介する.
\begin{equation*}
U^{n+1} = U^n + \Delta t \left( \frac{3}{2} f^n - \frac{1}{2} f^{n-1} \right).
\end{equation*}
右辺第 2 項の段階数を増やすことで精度を上げることができる. 打ち切り誤差は,
\begin{equation*}
f^{n-1} = f^n + \left( \DP{f}{U} \DP{U}{t} + \DP{f}{t} \right) \Delta t + O(\Delta t^2)
\end{equation*}
より,
\begin{align*}
\varepsilon_{AB} &= \frac{U^{n+1}-U^{n}}{\Delta t} - \left( \frac{3}{2} f^n - \frac{1}{2} f^{n-1} \right)\\
                 &= \left. \DD{U}{t} \right|_n + \frac{1}{2} \left. \DD[2]{U}{t} \right|_n \Delta t - \left\{ f^n - \frac{1}{2} \left( \DP{f}{U} \DP{U}{t} + \DP{f}{t} \right) \Delta t + O(\Delta t^2) \right\} + O(\Delta t^2)\\
                 &= \frac{1}{2} \DD{}{t} \left( \DD{U}{t} - f \right) \Delta t + O(\Delta t^2)\\
                 &= O(\Delta t^2).
\end{align*}

\section{振動方程式への応用}
\setcounter{equation}{0}
\begin{align}
\DD{U}{t} &= f(U,t),\\
f(U,t) &= i \omega U \ \ (\omega \in \mathbf{R}) \notag
\end{align}
について, 様々な時間差分スキームを適用し, その安定性を調べる. 様々な偏微分方程式は最終的に振動方程式を解く問題に帰着することが多い.

{\bf 例1) 一次元線形移流方程式}

\begin{equation*}
\DP{u}{t} + c \DP{u}{x} =0.
\end{equation*}
$\displaystyle u = Re \left[ U(t) e^{ikx} \right]$ とおくと,
\begin{equation*}
\DD{U}{t} + ikcU =0
\end{equation*}
となって, 振動方程式 (9.2.1) において $\omega = -kc$ とおいたものに等しくなる.

{\bf 例2) 慣性振動}

\begin{equation*}
\DD{u}{t} = fv \ \ , \ \ \DD{v}{t} = -fu.
\end{equation*}
複素速度 $U \equiv u+iv$ を導入すると,
\begin{align*}
\DD{u}{t} + i \DD{v}{t} &= fv -ifu\\
                        &= -if(u+iv).
\end{align*}
ゆえに,
\begin{equation*}
\DD{U}{t} = -ifU.
\end{equation*}
これは, $\omega = -f$ とした振動方程式である.

ちなみに, 振動方程式 (2.1) の一般解は,
\begin{equation*}
U(t) = U(0)e^{i\omega t},
\end{equation*}
であり, $t = n\Delta t$ の場合,
\begin{equation*}
U(n\Delta t) = U(0) e^{i\omega n \Delta t}
\end{equation*}
となる.

フォンノイマン法による安定性解析のため, 増幅係数を定義する.
\begin{equation}
\lambda \equiv \frac{U^{n+1}}{U^n}.
\end{equation}
ただし,
\begin{equation}
\lambda = \lvert \lambda \rvert e^{in\theta}.
\end{equation}
この時,
\begin{equation}
U^n = U^0 \lvert \lambda \rvert^n e^{in\theta}.
\end{equation}
安定性は以下の様に評価される.
\begin{alignat*}{2}
\lvert \lambda \rvert &> 1 & \qquad & \text{不安定}\\
\lvert \lambda \rvert &= 1 & \qquad & \text{中立}\\
\lvert \lambda \rvert &< 1 & \qquad & \text{減衰}
\end{alignat*}
位相については, 真の解と数値解との位相比をとって与える.
\begin{equation*}
\frac{in\theta}{i \omega n \Delta t} = \frac{\theta}{\omega \Delta t}.
\end{equation*}
評価は以下の通りである.
\begin{alignat*}{2}
\frac{\theta}{\omega \Delta t} &> 1 & \qquad & \text{位相は速く進む}\\
\frac{\theta}{\omega \Delta t} &= 1 & \qquad & \text{位相は一致}\\
\frac{\theta}{\omega \Delta t} &< 1 & \qquad & \text{位相は遅く進む}
\end{alignat*}
正確な数値解を得るには$\lvert \lambda \rvert $, $\displaystyle \frac{\theta}{\omega \Delta t}$ ともに 1 に近い方がよい. そうでない場合, 「計算モード」と呼ばれる偽りの解が現れることがある. この「計算モード」は$\Delta t , \Delta x \to 0$ にしても, 真の解に一致しない. 「計算モード」の振幅を抑制するためには, $\lvert \lambda \rvert < 1$ とした方がよい.

\subsection{反復しない 2 段階スキーム}
\begin{equation*}
U^{n+1} = U^n + \Delta t \left( \alpha f^n + \beta f^{n+1} \right).
\end{equation*}
但し,
\begin{equation*}
\alpha + \beta = 1.
\end{equation*}
である. つまり,
\begin{alignat*}{2}
\alpha &= 1 , \beta = 0 & \qquad & \text{オイラースキーム}\\
\alpha &= 0 , \beta = 1 & \qquad & \text{後退スキーム}\\
\alpha &= \beta = \frac{1}{2} & \qquad & \text{台形スキーム}
\end{alignat*}
振動方程式では, $f = i \omega U$ なので, 差分式は,
\begin{equation*}
U^{n+1} = U^n + i \omega \Delta t \left( \alpha U^n + \beta U^{n+1} \right).
\end{equation*}
$\displaystyle \lambda = \frac{U^{n+1}}{U^n}$ とすると,
\begin{equation*}
\lambda = 1 + i \omega \Delta t (\alpha + \beta \lambda ).
\end{equation*}
ゆえに,
\begin{align*}
\lambda &= \frac{1 + i\omega \Delta t \alpha}{1 - i\omega \Delta t \beta}\\
        &= \frac{1 + i\alpha p}{1 - i \beta p}.
\end{align*}
但し, $p \equiv \omega \Delta t$ である. よって,
\begin{equation*}
\lambda = \frac{1}{1+\beta^2 p^2} (1- \alpha \beta p^2 + ip).
\end{equation*}
$\alpha = 1, \beta = 0 $ (オイラースキーム) の場合,
\begin{equation*}
\lambda = 1+ip.
\end{equation*}
$\alpha = 0, \beta = 1 $ (後退スキーム) の場合,
\begin{equation*}
\lambda = \frac{1+ip}{1+p^2}.
\end{equation*}
$\alpha = \beta = \frac{1}{2} $ (台形スキーム) の場合,
\begin{equation*}
\lambda = \frac{1}{1 + \frac{1}{4} p^2} \left( 1-\frac{1}{4} p^2 + ip \right).
\end{equation*}
それぞれのスキームについて$\lvert \lambda \rvert$ を調べる.

{\bf オイラースキームの場合}

\begin{equation*}
\lvert \lambda \rvert = \sqrt{\lambda \lambda^*} = \sqrt{1+p^2} > 1.
\end{equation*}
よって, オイラースキームは振動方程式に対し不安定. 但し, $p = \omega \Delta t \ll 1$ のとき,
\begin{align*}
\lvert \lambda \rvert &= 1 + \frac{1}{2}p^2 - \frac{1}{4} p^4 + \cdots \\
                      &= 1 + O(\Delta t^2)\\
                      &= 1 + O(\Delta t^2),
\end{align*}
となり, フォンノイマン法の安定性の必要条件は満たしている.

{\bf 後退スキームの場合}

\begin{equation*}
\lvert \lambda \rvert = \sqrt{\lambda \lambda^*} = \frac{1}{1+p^2} \sqrt{1+p^2} = \frac{1}{1+p^2} < 1
\end{equation*}
よって, 後退スキームは $\Delta t$ の大きさによらず安定. 但し, $\omega $ が大きいほど減衰率も大きくなる. 実際の問題では, $\omega$ の大きい解 (高周波数解) は数値的に増幅しやすい (初期値の誤差のため). したがって, 後退差分スキームのような振動数によって選択的に減衰させるスキームは不要な高周波数解を除去するスキーム (フィルター) としても用いられる.

{\bf 台形スキームの場合}

\begin{align*}
\lvert \lambda \rvert &= \sqrt{\lambda \lambda^*} \\
                      &= \frac{1}{1+ \frac{1}{4} p^2} \sqrt{ \left(1 -\frac{1}{4} p \right)^2 + p}\\
                      &= \frac{1}{1+ \frac{1}{4} p^2} \sqrt{ \left(1 +\frac{1}{4} p \right)^2}\\
                      &= 1.
\end{align*}
よって, 台形スキームは中立である. 以上より, 陰的なスキームは$\Delta t$ の大きさによらず安定である.


\end{document}
