% 表題   DCPAM5  座標系・変換公式に関する解説 -- Gauss の公式, Gauss - Legendre函数の公式
%
% 履歴
%\Drireki{1991/12/04 保坂征宏}
%\Drireki{2005/04/04 石渡正樹}

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

  ここでは Gauss の台形公式を示す. 
  
  波数 $M$ 以下の三角函数で表現される 
  $g(\lambda)$ \ （$0 \le \lambda < 2 \pi$）
  \begin{align}
    g(\lambda) = \sum_{m=-M}^{m=M} g_m \exp( i m \lambda )
  \end{align}
  について $M < I$ を満たすように $I$ をとると, 
  \begin{align}
     &  \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) \nonumber
  \end{align}
  が成り立つ. これを Gauss の台形公式という. 

  より実用的な公式は, 
  \begin{align}
    &  \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) \nonumber
  \end{align}
  である. この証明は, 
  $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{align}
          \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{align}
   }. 

  
   以下に Gauss の台形公式の証明を記す. 
   まず, 左辺を計算すると, 
   \begin{align}
      \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{align}
   である. ここで, $\int_0^{2\pi} \exp(i m \lambda) d \lambda $ は
   $m=0$ の項しか残らないことを使った.
   一方右辺は
   \begin{align}
      \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) \nonumber \\
      &=  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{align}
   ここで, 上に示した「より実用的な公式」により
   \begin{align}
       \sum_{n=1}^{I} 
         \left( \exp(\frac{2 \pi i m}{I}) \right)^{n-1} =0
         \ \ \ \ (m \neq 0)
   \end{align}
   が成り立つ. したがって, 
   \begin{align}
      \frac{1}{2\pi} \int_0^{2\pi} g(\lambda) d \lambda 
        = \frac{1}{I} \sum_{n=1}^{I} g(\lambda_n) 
   \end{align}
   となる. 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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{align}
      &  \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{align}
   ここで, $w_j$ は Gauss 荷重と呼ばれる. 
   
   以下では上の式を証明する. 
   ただし, Legendre函数としては, 
   最初は岩波公式集のLegendre函数 $\tilde{P_n}$ を用い, 
   最後に2で規格化したLegendre函数 $P_n$ に
   直すことにする
          \footnote{
                混乱を招かぬよう, このような手続きを踏む. 
                実際, 公式集を含む他の文献には 
                $\tilde{P}_n^m$ の公式が書かれていることが多いので, 
                このように書く方が他と参照しやすいであろう. }. 
   

   \underline{STEP 1} \ \ Lagrange 補間の導入

   $f(\mu)$ を $K$ 次多項式（$0 \le K \le 2J-1$）とする. 
   $\tilde{P}_n$ を岩波公式集のLegendre函数（Rodriguesの公式）
   とする. 
   %
   \begin{align}
     \int_{-1}^1 \tilde{P}_n(\mu) \tilde{P}_{n'}(\mu) d \mu 
        = \frac{2}{2n+1} \delta_{nn'} .
   \end{align}
   %
   $L(\mu)$ を, 
   $f(\mu_j)$ を Lagrange 補間公式にしたがって補間した多項式として
   定義する. 
   %
   \begin{align}
     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{align}
   このとき, 各$j$ について $ L(\mu_j) = f(\mu_j) $ である. 
   ここで $L$ は, 
    $0 \le K \le J-1$ の時（$f$ が $J-1$ 次以下の多項式）のときは 
    厳密に $L=f$ になる
               \footnote{
                    このことは
                    $L-f$ が $J-1$ 次以下の多項式であること,
                    $J$ 個の零点 $\mu_j$ を持つことから明らか. }
   ことに注意せよ. 

   したがって, 
   関数 $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{align}
             f(\mu) - L(\mu) = \tilde{P}_J(\mu) S(\mu)
           \end{align}
           と書くことができる. 
   \end{itemize}
   
   $f(\mu)-L(\mu)$ を $\mu$ について $-1$ から 1 まで積分する. 
   $J \le K \le 2J-1$ の時については
   Legendre函数の直交性より, 
   $\tilde{P}_J(\mu) S(\mu)$ の積分は零である.
   したがって, 
   \begin{align}
     \int_{-1}^{1} f(\mu) d \mu
       &=  \int_{-1}^1 L(\mu) d \mu \nonumber \\
       &=  \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                            \nonumber    \\  
       &=  \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 \nonumber \\
       &=  2 \sum_{j=1}^{J} f(\mu_j) w_j   
   \end{align}
   
   ここで, 証明すべき式の $P_J$ は規格化されていて, 
   上の式の $\tilde{P}_J$ は規格化されていないのにもかかわらず
   同じ $w_j$ が使われているが, 
   $\tilde{P}_J$ と $P_J$ の規格化定数は同じなので
   consistent である. 
   
   \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{align}
     (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{align}
   より, 
   \begin{align}
     (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| 
     \nonumber \\
      &=  
          (2n+1)(x-y)\tilde{P}_n(x)\tilde{P}_n(y)        \nonumber  \\
      & \quad
              +n (- \tilde{P}_{n-1}(x) \tilde{P}_n(y) 
              + \tilde{P}_{n-1}(y) \tilde{P}_n(x) )      \nonumber  \\
      &=  
          (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{align}
   となる. 
   この式を $n=0,1,\cdots,n-1$ について加えると, 
   \begin{align}
          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{align}
   が成り立つ. 
   ここで $n=J,x=\mu,y=\mu_j$ とすると $\tilde{P}_J(\mu_j)=0$ より, 
   \begin{align}
     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).
   \end{align}
   よって, 
   \begin{align}
       \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{align}
   である. 
   したがって, 
   \begin{align}
     w_j &=  \frac{1}{2} \int_{-1}^1 
              \frac{\tilde{P}_J(\mu)} 
                   {(\mu-\mu_j) \tilde{P}^{'}_J(\mu_j)} 
              d \mu                                 \nonumber \\
         &=  \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 \Deqlab{w 積分表示}
                                                    \nonumber \\
         &=  \frac{1}{J \tilde{P}_{J-1}(\mu_j) 
                        \tilde{P}^{'}_{J} (\mu_j)}
   \end{align}
   である. 
   ただし, \Deqref{w 積分表示} における積分は, $k=0$ の時のみ0でない値を
   持つこと, および $\tilde{P}_0 = 1$ を使った.
   さらに, 漸化式
   \begin{align}
     (1-\mu^2) \DP{\tilde{P}_n}{\mu}
         = n \tilde{P}_{n-1}(\mu) - n \mu \tilde{P}_n(\mu) 
   \end{align}
   で $n=J, \mu=\mu_j$ とする. $\tilde{P}_J (\mu_j)=0$ より,  
   \begin{align}
     w_j &=  \frac{1-\mu_j^2}{(J \tilde{P}_{J-1}(\mu_j))^2 }
   \end{align}
   となる. 
   
   \underline{STEP3} \ $\tilde{P}_n$ の規格化
   
   $P_n$ を
   \begin{align}
     \int_{-1}^1 P_n(\mu) P_n'(\mu) d \mu = 2
   \end{align}
   になるように規格化する. 
   ${\displaystyle 
      \tilde{P}_{J-1}
      =\sqrt{ \frac{1}{2(\mbox{J}-1)+1} } P_{J-1} 
    }$ 
   より, 
   \begin{align}
     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{align}
   となる. 
   
   \underline{まとめ}
   
   以上より
   \begin{align}
      &  \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{align}