\documentclass[12pt]{jarticle}

\oddsidemargin 0mm
\evensidemargin 0mm
\textwidth 160mm
\topmargin -15mm
\textheight 240mm

\author{}
\title{数値計算法概論：No.6 数値積分(Richardson 加速)}
\date{}

\begin{document}

\maketitle

\section{数値積分の収束加速法}

数値積分で台形公式を使って、刻み幅を小さくして精度をあ
げようとすると計算量が増大する。又、無闇に刻み幅を小さくすると今度は丸
め誤差の影響がでてくる。そこで刻み幅の荒いデータから収束性を早くする方
法を考えよう。

\subsection{Euler-Maclaurin 総和公式}

$f(x)$ が滑らかな関数の時、次のオイラー・マクローリン( Euler-Maclaurin )
総和公式が成立する
\cite{Recipe,Sugiwara}。
\begin{eqnarray}
\int_{a}^{b} f(x) dx &=& h \left[\frac{1}{2} f(a) + \sum_{i=1}^{i=N-1} f(a+h i)
+\frac{1}{2} f(b) \right]	\nonumber\\
&-& \frac{B_2 h^2}{2!}(f'(b)-f'(a)) - \cdots 
-\frac{B_{2k} h^{2k}}{(2k)!}(f^{(2k-1)}(b)- f^{(2k-1)}(a)) - \cdots,	\nonumber\\
h &\equiv& \frac{(b-a)}{N}	
\label{eq:Euler-Maclaurin}
\end{eqnarray}
ここで $B_{2k}$ は Bernoulli 数と呼ばれ、次の関係で定義される。
\begin{equation}
	\frac{t}{e^t-1} = \sum_{n=0}^\infty B_n \frac{t^n}{n!}
\end{equation}
具体的には 
\begin{equation}
 B_0=1,B_2=\frac{1}{6}, B_4=-\frac{1}{30},B_6= \frac{1}{42}, \cdots
\end{equation}
(奇数値は $B_1=-1/2$ を除くとすべて消える)。

\subsection{収束性}

分割数 $N$ の時の台形公式による数値積分の結果を $I_N$ とする。
\begin{equation}
I_N = h_N \left[\frac{1}{2} f(a) + \sum_{i=1}^{i=N-1} f( a + h_N i)
+\frac{1}{2}f(b) \right] ; \;
h_{N} = \frac{b-a}{N}
\end{equation}
台形公式の収束性はオイラー・マクローリンの公式(\ref{eq:Euler-Maclaurin})
で与えられ、$h$ の偶関数で
ある。分割数 $N$ に対して誤差は $N^{-2},N^{-4},\ldots$ と振舞う。
\begin{equation}
 I_N= I_\infty + c_1 N^{-2} + c_2 N^{-4} + c_3 N^{-6} \cdots
\end{equation}

ところで $I_N$ と $I_\infty$ の直接比較して収束性評価できるのは結果が解
析的に分かっている場合のみである。
一般には、$I_{2N}-I_{N}$ を調べることで
収束性を調べられる。例えば誤差が、
\begin{equation}
 I_N = I_\infty + c N^{-\alpha}
\end{equation}
と冪乗的に振舞うなら、
$I_{2N}-I_{N}=c'N^{-\alpha},\quad c'=c (2^{-\alpha}-1)$ 
となる。
さらに
\begin{equation}
 (I_{4N}-I_{2N})/(I_{2N}-I_{N}) = 2^{-\alpha}
\end{equation}
として誤差の収束の早さも評価できる。

\subsection{台形公式の精度向上法:Richardson 補外}

さらにオイラー・マクローリンの公式
(\ref{eq:Euler-Maclaurin})
を使うと収束性を改善できる。つまり、
\[
 I_N= I_\infty + c_1 N^{-2} + c_2 N^{-4} + c_3 N^{-6} \cdots 
\]
と
\[
 I_{2N}= I_\infty + c_1 (2N)^{-2} + c_2 (2N)^{-4} + c_3 (2N)^{-6} \cdots 
\]
を比較して、
\begin{equation}
 I_N^{(1)} = \frac{I_{2N} -I_{N}/4}{1-1/4}
\end{equation}
とすると $N^{-2}$ の誤差の項が打ち消されるので、
\begin{equation}
 I_N^{(1)}- I_\infty = - \frac{c_2}{4} N^{-4}- \frac{5 c_3}{16} N^{-6}\cdots
\end{equation}
となり元の $I_N$ より収束が早くなる。この操作を Richardson 加速
\cite{Iri,Sugiwara2}と呼ぶ。台形公式による数
値積分を Richardson 加速すると、誤差は一般に $N^{-4}$ のように振舞う。

上のようにして得られた $I_N^{(1)}$ の数列を、さらに
\begin{equation}
 I_N^{(2)} = \frac{I_{2N}^{(1)} -I_{N}^{(1)}/16}{1-1/16}
\end{equation}
とすると、誤差が $N^{-6}$ のように振舞う。

特に分割数$N$を$2,4,8,16,32,\cdots$ と2の冪乗$N=2^m$にとると、
\begin{equation}
 I_{2^m} = I_\infty + \sum_{j=1}^\infty c_j 2^{- 2 j m}
\end{equation}
であるので、上の操作を繰り返すことができて、
\begin{equation}
 I_N^{(k)} = \frac{I_{2N}^{(k-1)} -I_{N}^{(k-1)}/4^{k}}{1-1/4^{k}}
\end{equation}
とするとさらに収束を早くできる。

\section{プログラム}

積分の上限、下限と分割数を与えた時、台形公式で数値積分の結果を出すサブルー
チンを利用すると、分割数を2の冪乗とした時の結果を出すプログラムは以下の
ようになる。
\begin{verbatim}
c     Main routine
c
      implicit real*8 (a-h,o-z)
      dimension dInt(10,3)
c
      do m=1,6
         n=2**m
         dInt(m,1)=dIntegral(n,a,b)
      end do
c
      do m=1,6-1
         write(*,*) 2**m, dInt(m,1),dInt(m+1,1)-dInt(m,1)
      end do
      end
c
c
c     Integral routine(trapezoidal integral)
c
      real*8 function dIntegral(n,a,b)
      implicit real*8 (a-h,o-z)
      integer n,i
c
      ...........
      dIntegral=......
      return
      end
\end{verbatim}
また、これを利用すると、Richardson 補外のプログラムができる。複数回 
Richardson 補外を繰り返すには、subroutine 文を使うと有効である。

\section{課題}

(1-1) No.5 の積分に付いて、分割数$N$を変えて $I_N$ を計算し
$O(h^2) (=O(1/N^2))$の誤差の収束性を確かめよ。
この中で、収束性が予期されるものよりずっと遅いものが
ある。その原因を考察して見よ。

(1-2) さらに Richardson 補外をして $I_N^{(1)},I_N^{(2)} \ldots$ の収束性
を見よ。
以上のものの内、収束性が予想より早いものと、ずっと遅いものと
がある。その原因を考察して見よ。

\begin{thebibliography}{99}
\bibitem{Recipe}「Numerical Recipes in C」、技術評論社、
第4章2節
\bibitem{Sugiwara}「数値計算法の数理」、杉原 正顕、室田 一雄 著、
	岩波書店、第12章3節
\bibitem{Iri}「数値計算の常識」、伊理 正夫、藤野 和建 著、
	共立出版、第7章
\bibitem{Sugiwara2}「数値計算法の数理」、杉原 正顕、室田 一雄 著、
	岩波書店、第12章4節
\end{thebibliography}

\end{document}

