\documentclass[a4j,12pt]{jreport}
\usepackage{Dennou6}
\usepackage{times}
\usepackage{ascmac}
\usepackage{float}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{bm}

\Dauthor{馬場 健聡}
\Ddate{2010 年 01 月 13 日}% 日付の指定
\Dtitle{GFDワークノート}      % 左上ヘッダー指定
\Dfile{2011-0113-kintore.tex}
\markright{第10章 1次元移流方程式の数値解法} % 右上ヘッダの指定

\pagenumbering{roman}
\setcounter{page}{6}       
\setcounter{chapter}{10}     
\setcounter{section}{1}     
\setcounter{subsection}{1}  
\setcounter{equation}{0}   
\setcounter{figure}{0}      
\setlength{\parskip}{1em}   

%\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{chapter}.\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}}
\def\thesubsection{\arabic{section}.\arabic{subsection}}

\begin{document}

\section{数値分散}
一次元線形移流方程式を数値的に解くと, 使用する空間差分法によっては「数値分散」
という現象が現れる. ここでは数値分散とその仕組みに関して解説する. 
元の方程式について考える. 
\begin{figure}[hb]
\begin{center}
\includegraphics[width=100mm]{./fig/fig1.eps}
\caption{
}
\end{center}
\end{figure}
\begin{align}
  \DP{u}{t} + c \DP{u}{x} = 0. \Deqlab{1} 
\end{align}
$u(x,t)$ をフーリエ展開し, そのとある波数成分の解, 
\begin{align}
  u(x,t) = Re [U(t) e^{ikx}] \Deqlab{2}
\end{align}
について考える. これを\Deqref{1}に代入すると
$U(t)$ の満たす式
\begin{align}
  \DD{U(t)}{t} + ikc U(t) = 0 \Deqlab{3}
\end{align}
を得る.
この式は振動方程式
\begin{align}
  \DD{U(t)}{t} + \nu U(t) = 0 \nonumber
\end{align}
と同じ形をしている.　 
位相速度 $c_{ph}$ は定義から, 
\begin{align}
c_{ph} \equiv \frac{\nu}{k} = \frac{kc}{k} = c. \nonumber
\end{align}

以上のように, 元の式 (連続形の式)では分散性はない. 
(位相速度は波形によらない.) \\

次に, 空間微分のみ差分化した式
\begin{align}
\DP{u_j(t)}{t} + c \DP{u_{j+1} - u_{j-1}}{2 \Delta x} = 0 \Deqlab{4}
\end{align}
について考える. ここで, 空間微分を中心差分で表した\Deqref{4}の
時間微分を適当な方法で差分化すると, 元の\Deqref{1}の有限差分法
が得られる. 

ex)オイラー法を用いると\\
\begin{align}
 \frac{u_j^{n+1} - u_j^n}{\Delta t} + \frac{u_{j+1}^n - u_{j-1}^n}{\Delta x} = 0 \Deqlab{5}
\end{align}
\Deqref{4}は\Deqref{5}における $\Delta t \rightarrow 0$ の極限になっている. \\

連続形の式\Deqref{1}について考察したのと同様に, \Deqref{4}に対し, 
\begin{align}
u_j(t) = Re[U_j(t) e^{ikj \Delta x}] \nonumber
\end{align}
を代入する. この時, $U_j(t)$ の満たす式として
\begin{align}
\DD{U_j}{t} + ik \left( c \frac{\sin k \Delta x}{k \Delta x} \right) u_j =0 \Deqlab{6}
\end{align}
$\nu^* = k \left( c \frac{\sin k \Delta x}{k \Delta x} \right)$
 と置くと位相速度 $c^*_ph$ は
\begin{align}
c^*_ph = \frac{\nu^*}{k} = c \frac{\sin k \Delta x}{k \Delta x} \Deqlab{7}
\end{align}
となる. 位相速度は波数 $k$ に依存するため, \Deqref{4}の解には分散性がある. 
離散化したことにより, 解に分散性が現れることを数値分散とよぶ. \\
\Deqref{4}の $k$ 依存性を調べる. 
$k\Delta x$ が 0 から増加するにつれ $c_{ph}^*$ は単調減少する. 
解像可能な最小波長の成分 $(\lambda = 2 \Delta x)$ に対し, 
$c_{ph}^* =0$ になる. ($k=\frac{2\pi}{\lambda}=\frac{\pi}{\Delta x}$)\\

まとめると \\
・どの波数成分の位相速度も解析値よりも小さい \\
・$\lambda = 2\Delta x$ の成分の位相速度は 0 (定常波になる) \\

\begin{figure}[hbtp]
\begin{center}
\includegraphics[width=100mm]{./fig/fig2.eps}
\caption{
}
\end{center}
\end{figure}

次に\Deqref{7}の解の群速度を計算する. 
定義から群速度 $c_g$ は
\begin{align}
c_g = \DD{\nu}{k} \nonumber
\end{align}
連続形の式\Deqref{3}では $\nu = kc$ なので, 
\begin{align}
c_g \equiv  \DD{kc}{k} = c \nonumber
\end{align}
となる. 
一方, 
\begin{align}
  \nu^* &= kc^* = k \left( c \frac{\sin k \Delta x}{ k \Delta x} \right) \nonumber
\end{align}
より, \Deqref{6}の解の群速度 $c_g^*$ は
\begin{align}
  c_g^* &= \DD{}{k} \left( c \frac{\sin k \Delta x}{ k \Delta x} \right) \nonumber \\
        &= c \cos k\Delta x \Deqlab{9}
\end{align}
となる. 

$k\Delta x$ が　0 から増加するとともに $c_g^*$ は単調減少し, \\
$k\Delta x = \frac{\pi}{2} (k=4 \Delta x)$ の時, $c_g =0$ \\
$k\Delta x = \pi (k=2 \Delta x)$ の時, $c_g =0$ となる. 

\subsubsection*{群速度の分散性の具体例}
空間的に滑らかな関数 $Y(x)$ を考える. (例えば波長の長い正弦波)
$Y(x)$ を空間方向に離散化したものを $Y_j$ とし, 
\begin{align}
Y_j = Y(x) \nonumber
\end{align}
とおく. 
\begin{figure}[hbtp]
\begin{center}
\includegraphics[width=100mm]{./fig/fig3.eps}
\caption{
}
\end{center}
\end{figure}\\

次に $Y_j$ を用いて表される
\begin{align}
u_j = (-1)^j Y_j \Deqlab{10}
\end{align}
という波を考える. これを\Deqref{4}式
\begin{align}
\DP{u_j}{t} + c \frac{u_{j+1} - u_{j-1}}{2 \Delta x} = 0 \nonumber
\end{align}
へ代入すると, 
\begin{align}
(-1)^j \DP{Y_j}{t} + c \frac{(-1)^{j+1} Y_{j+1} - (-1)^{j-1} Y_{j-1}}{2 \Delta x} &= 0 \nonumber \\
\DP{Y_j}{t} - c \frac{Y_{j+1} - Y_{j-1}}{2 \Delta x} &= 0 \nonumber
\end{align}
したがって, \Deqref{10} で定義された $u_j$ は $-c$ で移流される. \\
($x$ の負の方向に移流される.)\\

\Deqref{4} に戻り, その解析解を考える. 
  \begin{align}
    \DP{u_j}{t} + c\frac{u_{j+1} - u_{j-1}}{2 \Delta x} = 0 \nonumber
  \end{align}
無次元座標 $\tau$ を導入する. 
  \begin{align}
    \tau = \frac{ct}{\Delta x} \Deqlab{11}
  \end{align}
\Deqref{4}を $\frac{c}{2\Delta t}$ で割り, 
  \begin{align}
    \DP{u_j}{\left( \frac{c}{2\Delta x} \right)} + (u_{j+1} - u_{j-1}) &= 0 \nonumber \\
    2 \DP{u_j}{\tau} &= u_{j-1}(\tau) - u_{j+1}(\tau) \Deqlab{12}
  \end{align}
\Deqref{12}は第一種ベッセル関数の満たす漸化式に等しい. 
j 次の第一種ベッセル関数を $J_j$ と表すと, 
\begin{align}
u_j (\tau) = J_j(\tau ) \Deqlab{13}
\end{align}
\Deqref{4}において $j=0$ とする点は任意にとることが出来るので一般的には
\begin{align}
u_j (\tau) = J_{j-p}(\tau ) \nonumber
\end{align}
と表される. ここで $p$ は任意の整数である.
\Deqref{4}は線形方程式なので一般解は, 
\begin{align}
  u_j (\tau) = \sum_{p= - \infty}^{\infty} a_p J_{j-p} (\tau). \Deqlab{14}
\end{align}
$\tau=0$ のとき, $J_0$ を除く全ての $J_k$ の値は定義より $0$ になる. 
$J_0(\tau=0) = 1$ より, $\tau=0$ を\Deqref{14}に代入して, 
\begin{align}
u_j(0) = a_j.
\end{align}
$a_j$ は初期条件 $u_j(0)$ より与えられる. 
具体例として次のような初期条件を与えた場合を考える. 
\begin{eqnarray}
u_j(0) = \left\{
\begin{array}{l}
1 (j=0) \nonumber \\
0 (j \neq 0) \nonumber 
\end{array}
\right. \nonumber
\end{eqnarray} 

\begin{figure}[hbtp]
\begin{center}
\includegraphics[width=100mm]{./fig/fig4.eps}
\caption{一番上の図を初期条件とした時の\Deqref{4}の解析解. 
}
\end{center}
\end{figure}

初期条件 $u(x)$ が $\delta$ 関数の場合, 
これをフーリエ余弦展開すると, 
\begin{align}
u(x) &= \int_0^\infty a(k) \cos kx dk \nonumber \\
a(k) &= \frac{2}{\pi} \int_0^\infty u(x) \cos kx dx \nonumber
\end{align}
$\delta$ 関数は
\begin{eqnarray}
\delta (x) = \left\{
\begin{array}{l}
1 (x=0) \nonumber \\
0 (x \neq 0) \nonumber 
\end{array}
\right. \nonumber
\end{eqnarray} 
より $a(k)$ を数値的に求めると, 

\begin{align}
  a(k) &= \frac{2}{\pi} u(j\Delta x) \cos k(j\Delta x) \Delta x \nonumber \\
       &= \frac{2}{\pi} u(0) \cos (0) \Delta x \nonumber \\
       &= \frac{2}{\pi} \Delta x \nonumber
\end{align}

光学のアナロジーからこの様な擾乱は白色雑音 (white noize) と呼ばれる.

\end{document}
