% 表題   DCPAM2 第2部 離散化 Gauss の公式, Gauss - Legendre函数の公式
%
% 履歴
\Drireki{91/12/04 保坂征宏}
\Drireki{2005/04/04 石渡正樹}

\section{積分評価}
   
   
\subsection{Gauss の台形公式}

  ここでは Gauss の台形公式を示す. 
  \vspace{5mm}
  
  波数 $M$ 以下の三角函数で表現される 
  $g(\lambda)$ \ （$0 \le \lambda < 2 \pi$）
  \begin{eqnarray}
    g(\lambda) = \sum_{m=-M}^{m=M} g_m \exp( i m \lambda )
  \end{eqnarray}
  について $M < I$ を満たすように $I$ をとると, 
  \begin{eqnarray}
     & &\frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda 
       = \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n) \\
     & & \lambda_n = \frac{2\pi (n-1)}{I} \ (n=1,2,\cdots,I)
  \end{eqnarray}
  が成り立つ. これを Gauss の台形公式という. 
  \vspace{5mm}
  
  より実用的な公式は, 
  \begin{eqnarray}
    & &\sum_{n=1}^{I} \exp(i m \lambda_n) 
          = \left\{
                \begin{array}{ll}
                    I & \ \ m=0 の時 \\
                    0 & \ \ 0 < |m| < I の時 \\
                \end{array}
            \right.                                     \\
     & & \lambda_n = \frac{2\pi (n-1)}{I} \ (n=1,2,\cdots,I)
  \end{eqnarray}
  である. この証明は, 
  $I > M (|m| の最大値)$ より $m \neq 0$ の時には
  ${\displaystyle \exp(i m \lambda_n) 
                     = \exp \left( \frac{2\pi i m (n-1)}{I} \right) }$ 
  において, 全ての $n$ について 
  $m(n-1)$ が$I$ の整数倍になることがないことを考慮すると明らかである
  ( $m$, $n-1$ はともに $I$ よりも小さい整数なので, $m(n-1)$ は 
    $I$ の整数倍にならない)
  \footnote{
      等比級数の和を直接計算しても良い.
      \begin{eqnarray}
          \sum^{I}_{n=1} \exp \left\{ im \frac{2 \pi (n-1)}{I} \right\}
         = \frac{1 - \left( e^{\frac{im 2 \pi}{I}}\right)^I}
                {1 - e^{ \frac{im 2 \pi}{I} } }
         = \frac{1 -  e^{im 2 \pi} }
                {1 - e^{ \frac{im 2 \pi}{I} } }
         = 0
      \end{eqnarray}
   }. 
  \vspace{5mm}
  
  以下に Gauss の台形公式の証明を記す. \\
   まず, 左辺を計算すると, 
   \begin{eqnarray}
      \frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda 
         &=& \sum_{m=-M}^{M} \frac{1}{2\pi} g_m 
                \int_0^{2\pi} \exp(i m \lambda) d \lambda 
          =g_0 
   \end{eqnarray}
   である. ここで, $\int_0^{2\pi} \exp(i m \lambda) d \lambda $ は
   $m=0$ の項しか残らないことを使った.
   一方右辺は
   \begin{eqnarray}
      \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n) 
      &=& \frac{1}{I} 
           \sum_{n=1}^{I} \sum_{m=-M}^{M} 
           g_m \exp(i m \lambda_n) \\
      &=& g_0 
          + \sum_{m=-M, m\neq0}^{M} \frac{g_m}{I} 
             \sum_{n=1}^{I} 
               \left( \exp(\frac{2 \pi i m}{I}) \right)^{n-1} 
   \end{eqnarray}
   ここで, 上に示した「より実用的な公式」により
   \begin{eqnarray}
       \sum_{n=1}^{I} 
         \left( \exp(\frac{2 \pi i m}{I}) \right)^{n-1} =0
         \ \ \ \ (m \neq 0 \ の時)
   \end{eqnarray}
   が成り立つ. したがって, 
   \begin{eqnarray}
      \frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda 
        = \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n) 
   \end{eqnarray}
   となる. 
  


\subsection{Gauss-Legendreの公式}

   $f(\mu)$ を $2J-1$ 次以下の多項式とする.  
   $P_n$ を2で規格化した n 次の Legendre函数とする. 
   このとき, ${\displaystyle \int_{-1}^1 f d \mu }$ は 
   $P_J$ の零点である Gauss 格子 $\mu_j（j=1,2,\cdots,J）$ における 
   $f$ の値 $f(\mu_j)$ のみを用いて, 
   次式にもとづいて正確に評価することができる. 
   \begin{eqnarray}
      & &\int^{1}_{-1} f(\mu) d \mu = 2 \sum_{j=1}^{J} f(\mu_j) w_j  \\
      & &      w_j = \frac{1}{2} \int_{-1}^{1} 
                      \frac{P_J(\mu)}{(\mu-\mu_j) P^{'}_J (\mu_i)}d \mu
                   = \frac{(2J-1)(1-\mu_j^2)}{(J P_{J-1}(\mu_j))^2 }  .
   \end{eqnarray}
   ここで, $w_j$ は Gauss 荷重と呼ばれる. 
   \vspace{10mm}
   
   以下では上の式を証明する. 
   ただし, Legendre函数としては, 
   最初は岩波公式集のLegendre函数 $\tilde{P_n}$ を用い, 
   最後に2で規格化したLegendre函数 $P_n$ に
   直すことにする\footnotemark. 
          \footnotetext{
                私が混乱しないためにこのような手続きを踏む. 
                実際, 公式集を含む他の文献には 
                $\tilde{P}_n^m$ の公式が書かれていることが多いので, 
                このように書く方が他と参照しやすいであろう. }
   \vspace{5mm}
   
   \underline{STEP 1} \ \ Lagrange 補間の導入

   $f(\mu)$ を $K$ 次多項式（$0 \le K \le 2J-1$）とする. 
   $\tilde{P}_n$ を岩波公式集のLegendre函数（Rodriguesの公式）
   とする. 
   
   \begin{eqnarray}
     \int_{-1}^1 \tilde{P}_n(\mu) \tilde{P}_{n'}(\mu) d \mu 
        = \frac{2}{2n+1} \delta_{nn'}
   \end{eqnarray}
   
   $L(\mu)$ を, 
   $f(\mu_j)$ を Lagrange 補間公式にしたがって補間した多項式として
   定義する. 
   
   \begin{eqnarray}
     L(\mu) \equiv 
              \sum_{j=1}^{J} f(\mu_j) 
              \prod_{k=1,k \neq j}^{J} 
               \frac{ \mu-\mu_k}{\mu_j-\mu_k} 
   \end{eqnarray}
   このとき, 各$j$ について $ L(\mu_j) = f(\mu_j) $ である. 
   ここで $L$ は, 
    $0 \le K \le J-1$ の時（$f$ が $J-1$ 次以下の多項式）のときは 
    厳密に $L=f$ になる\footnotemark 
               \footnotetext{
                    このことは
                    $L-f$ が $J-1$ 次以下の多項式であること,
                    $J$ 個の零点 $\mu_j$ を持つことから明らか. }
   ことに注意せよ. 
   \vspace{5mm}
   
   したがって, 
   関数 $f(\mu)-L(\mu)$ は
   \begin{itemize}
     \item $0 \le K \le J-1$ の時, $0$ である. 
     \item $J \le K \le 2J-1$ の時, \\
           $\mu=\mu_j$ を零点とする $K$ 次多項式である. 
           $\mu_j$ は $J$ 次多項式 $\tilde{P}_J(\mu)$ の
           零点であることを思い出すと, 
           $f-L$ は $\tilde{P}_J(\mu)$ で割り切れるので, 
           ある $K-J$ 次多項式 $S(\mu)$ を用いて, 
           \begin{eqnarray}
             f(\mu) - L(\mu) = \tilde{P}_J(\mu) S(\mu)
           \end{eqnarray}
           と書くことができる. 
   \end{itemize}
   
   $f(\mu)-L(\mu)$ を $\mu$ について $-1$ から 1 まで積分する. 
   $J \le K \le 2J-1$ の時については
   Legendre函数の直交性より, 
   $\tilde{P}_J(\mu) S(\mu)$ の積分は零である.
   したがって, 
   \begin{eqnarray}
     \int_{-1}^{1} f(\mu) d \mu
       &=& \int_{-1}^1 L(\mu) d \mu \\
       &=& \sum_{j=1}^{J} f(\mu_j) 
            \int^1_{-1}
             \frac{ {\displaystyle
                        \prod_{k=1}^{J}(\mu-\mu_k) }
                  }
                  { (\mu-\mu_j)
                    {\displaystyle
                        \prod_{k=1,k \neq j}^{J}(\mu_j-\mu_k) }
                  }    
             d \mu                                         \\  
       &=& \sum_{j=1}^{J} f(\mu_j) 
             \int_{-1}^1
             \frac{\tilde{P}_J(\mu)}
                  {(\mu-\mu_j)\tilde{P}^{'}_J(\mu_j)} d \mu \\
       &=& 2 \sum_{j=1}^{J} f(\mu_j) w_j   
   \end{eqnarray}
   
   ここで, 証明すべき式の $P_J$ は規格化されていて, 
   上の式の $\tilde{P}_J$ は規格化されていないのにもかかわらず
   同じ $w_j$ が使われているが, 
   $\tilde{P}_J$ と $P_J$ の規格化定数は同じなので
   consistent である. 
   \vspace{5mm}
   
   \underline{STEP 2} \ \ 
   ${\displaystyle w_j
         = \frac{1}{2} \int_{-1}^1
             \frac{\tilde{P}_J(\mu)}
                  {(\mu-\mu_j)\tilde{P}^{'}_J(\mu_j)} d \mu
   }$ の漸化式を用いた変形
   
   漸化式 (岩波の Lgendre 関数・陪関数の従う漸化式) において
   $m=0$ とした式
   \begin{eqnarray}
     (n+1) \tilde{P}_{n+1}(\mu) 
       = (2n+1) \mu \tilde{P}_n(\mu) - n \tilde{P}_{n-1}(\mu)
       \ \ \ \ (n=0,1,2,\cdots)
   \end{eqnarray}
   より, 
   \begin{eqnarray}
     (n+1) \left| 
             \begin{array}{ll}
                \tilde{P}_{n+1}(x) & \tilde{P}_n(x) \\
                \tilde{P}_{n+1}(y) & \tilde{P}_n(y) 
             \end{array}
           \right| 
      &=& 
          \left| 
             \begin{array}{ll}
                (2n+1)x \tilde{P}_n(x)-n\tilde{P}_{n-1}(x) 
                   &   \tilde{P}_n(x) \\
                (2n+1)y \tilde{P}_n(y)-n\tilde{P}_{n-1}(y) 
                   &   \tilde{P}_n(y) 
             \end{array}
           \right| 
    \\
      &=& 
          (2n+1)(x-y)\tilde{P}_n(x)\tilde{P}_n(y)                           \\
          & & +n (- \tilde{P}_{n-1}(x) \tilde{P}_n(y) 
              + \tilde{P}_{n-1}(y) \tilde{P}_n(x) ) \\
      &=& 
          (2n+1)(x-y)\tilde{P}_n(x)\tilde{P}_n(y) 
          +
           n \left| 
                   \begin{array}{ll}
                      \tilde{P}_{n}(x) & \tilde{P}_{n-1}(x) \\
                      \tilde{P}_{n}(y) & \tilde{P}_{n-1}(y) 
                   \end{array}
                 \right| 
   \end{eqnarray}
   となる. 
   この式を $n=0,1,\cdots,n-1$ について加えると, 
   \begin{eqnarray}
          n \left| 
             \begin{array}{ll}
                \tilde{P}_n(x) & \tilde{P}_{n-1}(x) \\
                \tilde{P}_n(y) & \tilde{P}_{n-1}(y) 
             \end{array}
           \right| 
      &=&  \sum_{k=0}^{n-1} (2k+1)(x-y)\tilde{P}_k(x)\tilde{P}_k(y) 
   \end{eqnarray}
   が成り立つ. 
   ここで $n=J,x=\mu,y=\mu_j$ とすると $\tilde{P}_J(\mu_j)=0$ より, 
   \begin{eqnarray}
     J  \tilde{P}_J (\mu)  \tilde{P}_{J-1} (\mu_j) 
       &=&  \sum_{k=0}^{J-1} (2k+1) (\mu-\mu_j) 
              \tilde{P}_k(\mu) \tilde{P}_k(\mu_j) \\
     ∴  \ \ 
       \frac{\tilde{P}_J (\mu) }{\mu-\mu_j}
       &=& \frac{ 
                   {\displaystyle
                      \sum_{k=0}^{J-1} (2k+1) 
                      \tilde{P}_k(\mu) \tilde{P}_k(\mu_j) }
                 }
                 {J \tilde{P}_{J-1}(\mu_j)} 
   \end{eqnarray}
   である. 
   したがって, 
   \begin{eqnarray}
     w_j &=& \frac{1}{2} \int_{-1}^1 
              \frac{\tilde{P}_J(\mu)} 
                   {(\mu-\mu_j) \tilde{P}^{'}_J(\mu_j)} 
              d \mu                                            \\
         &=& \frac{1}{2 J \tilde{P}_{J-1}(\mu_j) 
                        \tilde{P}^{'}_{J} (\mu_j)}
             \sum_{k=0}^{J-1} 
                  (2k+1) \tilde{P}_k(\mu_j) 
                  \int_{-1}^1 \tilde{P}_k(\mu) d \mu \\
         &=& \frac{1}{J \tilde{P}_{J-1}(\mu_j) 
                        \tilde{P}^{'}_{J} (\mu_j)}      \\
         & &  （前式の積分は, k=0 の時のみ0でない値を持ち, 
                \ \tilde{P}_0 = 1.）
   \end{eqnarray}
   である. 
   さらに, 漸化式
   \begin{eqnarray}
     (1-\mu^2) \DP{\tilde{P}_n}{\mu}
         = n \tilde{P}_{n-1}(\mu) - n \mu \tilde{P}_n(\mu) 
   \end{eqnarray}
   で $n=J, \mu=\mu_j$ とする. $\tilde{P}_J (\mu_j)=0$ より,  
   \begin{eqnarray}
     w_j &=& \frac{1-\mu_j^2}{(J \tilde{P}_{J-1}(\mu_j))^2 }
   \end{eqnarray}
   となる. 
   \vspace{5mm}
   
   \underline{STEP3} \ $\tilde{P}_n$ の規格化
   
   $P_n$ を
   \begin{eqnarray}
     \int_{-1}^1 P_n(\mu) P_n'(\mu) d \mu = 2
   \end{eqnarray}
   になるように規格化する. 
   ${\displaystyle 
      \tilde{P}_{J-1}
      =\sqrt{ \frac{1}{2(\mbox{J}-1)+1} } P_{J-1} 
    }$ 
   より, 
   \begin{eqnarray}
     w_j = \frac{1-\mu_j^2}
                  {(J  \sqrt{ \frac{1}{2J-1} } P_{J-1}(\mu_j))^2 } 
         = \frac{(2J-1)(1-\mu_j^2)}
                  {(J  P_{J-1}(\mu_j))^2 }
   \end{eqnarray}
   となる. 
   \vspace{5mm}
   
   \underline{まとめ}
   
   以上より
   \begin{eqnarray}
      & &\int^{1}_{-1} f(\mu) d \mu = 2 \sum_{j=1}^{J} f(\mu_j) w_j  \\
      & &w_j = \frac{(2J-1)(1-\mu_j^2)}{(J P_{J-1}(\mu_j))^2 }
   \end{eqnarray}