台形公式（1次多項式），Newton-Cotesの式（Lagrange多項式），Gauss-Legendre求積

# 初めに
次のような求積を考えている．
```math
\begin{align}
    \int_{a}^{b} 
        f (x)
    \; dx
    &\approx
    \sum_{i=1}^{n} 
        w_i
        f (x_i)
\end{align}
```
参考: 

* [いまさら聞けない計算力学の定石](https://www.maruzen-publishing.co.jp/book/b10112229.html)
* [ラグランジュの補間公式とその応用例](https://manabitimes.jp/math/726)
* [ガウス求積法（Gauss–Legendre 公式）](https://manabitimes.jp/math/2671)
* [数値積分法](https://slpr.sakura.ne.jp/sikinote/docs/numeric/integration/)
* [FastGaussQuadrature.jlで数値積分しましょう](https://zenn.dev/hyrodium/articles/a7af62af417cc6)

# Methods
## 台形公式（1次多項式）
次の方法
```math
\begin{align}
    \int_{a}^{b} 
        f (x)
    \; dx
    &\approx
    \frac{b - a}{2} (f (a) + f (b))
    =: I_{\text{Trapz}} [f; a, b]
\end{align}
```
は，$f (x)$が1次以下の多項式なら等号で結ぶことができる．これを，台形近似と呼ぶ．

## Newton-Cotesの式（Lagrange多項式）
相異なる$n$個だけの点$\lbrace x_i, f (x_i) \rbrace_{i=1}^{n}$が与えられたとき，
```math
\begin{align}
    L_n (x_i)
    &= f (x_i)
\end{align}
```
を満たす$(n-1)$次以下の多項式$L_{n-1}$が一つ定まり，
```math
\begin{align}
    L_{n-1} (x)
    &= \sum_{i=1}^{n} f (x_i) l_i (x)
\end{align}
```
で与えられる．これをLagrange多項式と呼ぶ．$l_i (x)$は
```math
\begin{align}
    l_i (x)
    &= \prod_{j = 1; j \neq i}^{n} \frac{x   - x_j}{x_i - x_j}
\end{align}
```
である．ところで，Lagrange多項式を再掲する．
```math
\begin{align}
    L_{n-1} (x)
    &= \sum_{i=1}^{n} f (x_i) l_i (x)
\end{align}
```
これを，区間$[a, b]$で積分する．
```math
\begin{align}
    \int_{a}^{b}
        L_{n-1} (x)
    \; dx
    &= \int_{a}^{b} \sum_{i=1}^{n}
        f (x_i) l_i (x)
    \; dx
    \\
    &= \sum_{i=1}^{n} \int_{a}^{b}
        l_i (x)
    \; dx
    \; f (x_i) 
    \\
    &= \sum_{i=1}^{n} 
        w_i f (x_i) 
    \quad
    \left(
        w_i = \int_{a}^{b} l_i (x) \; dx
    \right)
\end{align}
```
$L_{n-1}$は$(n-1)$次以下の多項式なのだから，
```math
\begin{align}
    \int_{a}^{b}
        f (x)
    \; dx
    &\approx
    \sum_{i=1}^{n} 
        w_i f (x_i) 
    =: I_{\text{NC}} [f; a, b, n]
    \quad
    \left(
        w_i = \int_{a}^{b} l_i (x) \; dx
    \right)
\end{align}
```
は，$f (x)$が$(n-1)$次以下の多項式なら等号で結ばれる．これを，Newton-Cotesの式と呼ぶ（特に$x_i$を等幅に選んだものをNewton-Cotesの式と呼ぶらしいが，何れにせよ，$(n-1)$次以下の多項式で厳密という性質に変わりはない）．

## Gauss-Legendre求積
以降，積分区間が$[-1, 1]$である場合を考える．特に，$x \in [a, b]$である場合には，
```math
\begin{alignat}{3}
    x
    &= \frac{b - a}{2} y + \frac{b + a}{2}
    \quad
    &&\Leftrightarrow
    \quad
    y
    &&= \frac{2}{b - a}(x - a) - 1
    \\
    dx
    &= \frac{b - a}{2} dy
    \quad
    &&\Leftrightarrow
    \quad
    dy
    &&= \frac{2}{b - a} dx
\end{alignat}
```
とすることで，$y \in [-1, 1]$を考えることができる．このとき，$n$次のLegendre多項式$P_{n}$の零点を積分点に採用した
```math
\begin{align}
    \int_{-1}^{1}
        f (x)
    \; dx
    &\approx
    \sum_{i=1}^{n} 
        w_i f (x_i) 
    =: I_{\text{GL}} [f; n]
    \quad
    \left(
        w_i = \int_{-1}^{1} l_i (x) \; dx,
        \;
        P_n (x_i) = 0
    \right)
\end{align}
```
は，$f (x)$が$(2 n-1)$次以下の多項式なら等号で結ばれる．これを，Gauss-Legendre求積と呼ぶ．

# Results
## 3次多項式
[Wikipedia](https://en.wikipedia.org/wiki/Gaussian_quadrature)に記載の，
```math
\begin{align}
    f (x)
    &= 7 x^3 - 8 x^2 - 3 x + 3
\end{align}
```
を$[-1, 1]$で積分する．厳密解は$2 / 3$．

**台形則**
端点を集める．
```math
\begin{align}
    \int_{a}^{b} 
        f (x)
    \; dx
    \approx
    I_{\text{Trapz}} [f]
    &=
    \frac{b - a}{2} (f (a) + f (b))
    \\
    &= \frac{1 - (-1)}{2} (f (-1) + f (1))
    \\
    &= - 10
\end{align}
```
poor過ぎる．

**Newton-Cotes**
再掲すると，
```math
\begin{align}
    \int_{a}^{b}
        f (x)
    \; dx
    \approx
    I_{\text{NC}} [f; n]
    &=
    \sum_{i=1}^{n} 
        w_i f (x_i) 
    \\
    w_i
    &= \int_{a}^{b} l_i (x) \; dx
    \\
    l_i (x)
    &= \prod_{j = 1; j \neq i}^{n} \frac{x   - x_j}{x_i - x_j}
\end{align}
```
である．

$n = 2$として，$(x_1, x_2) = (-0.5, 0.5)$を採用してみる（端点を選べば台形公式に一致）．
```math
\begin{align}
    w_1
    &= \int_{-1}^{1} l_1 (x) \; dx
     = \int_{-1}^{1} \left( \frac{1}{2} - x \right) \; dx
     = 1
    \\
    w_2
    &= \int_{-1}^{1} l_2 (x) \; dx
     = \int_{-1}^{1} \left( \frac{1}{2} + x \right) \; dx
     = 1
    \\
    I_{\text{NC}} [f; n=2]
    &= \sum_{i=1}^{n=2} w_i f (x_i) 
     = - 2
\end{align}
```
未だ非常にpoor．

$n = 3$にしてみよう（$(x_1, x_2, x_3) = (-0.5, 0, 0.5)$）．
```math
\begin{align}
    w_1
    &= \int_{-1}^{1} l_1 (x) \; dx
     = \frac{4}{3}
    \\
    w_2
    &= \int_{-1}^{1} l_2 (x) \; dx
     = - \frac{2}{3}
    \\
    w_3
    &= \int_{-1}^{1} l_3 (x) \; dx
     = \frac{4}{3}
    \\
    I_{\text{NC}} [f; n=3]
    &= \sum_{i=1}^{n=3} w_i f (x_i) 
     = \frac{2}{3}
\end{align}
```
正確になった．いまの多項式$f = 7 x^3 - 8 x^2 - 3 x + 3$は3次であるが，$x^3$を積分すると落ちている．

$n = 3$のまま，$(x_1, x_2, x_3) = (-1, 0, 1)$に改めてみる．
```math
\begin{align}
    w_1
    &= \int_{-1}^{1} l_1 (x) \; dx
     = \frac{1}{3}
    \\
    w_2
    &= \int_{-1}^{1} l_2 (x) \; dx
     = \frac{4}{3}
    \\
    w_3
    &= \int_{-1}^{1} l_3 (x) \; dx
     = \frac{1}{3}
    \\
    I_{\text{NC}} [f; n=3]
    &= \sum_{i=1}^{n=3} w_i f (x_i) 
     = \frac{2}{3}
\end{align}
```
$n=3$個だけの積分点があれば，結果に影響はしないようである．Lagrange多項式を考えれば当然だろうか．

**Gauss-Legendre**
再掲である．
```math
\begin{align}
    \int_{-1}^{1}
        f (x)
    \; dx
    \approx
    I_{\text{GL}} [f; n]
    &=
    \sum_{i=1}^{n} 
        w_i f (x_i) 
    \\
    w_i
    &= \int_{-1}^{1} l_i (x) \; dx
    \\
    l_i (x)
    &= \prod_{j = 1; j \neq i}^{n} \frac{x   - x_j}{x_i - x_j}
    \\
    P_n (x_i)
    &= 0
\end{align}
```
$n = 2$とする．
```math
\begin{align}
    w_1
    &= \int_{-1}^{1} l_1 (x) \; dx
     = 1
    \\
    w_2
    &= \int_{-1}^{1} l_2 (x) \; dx
     = 1
    \\
    I_{\text{GL}} [f; n=2]
    &= \sum_{i=1}^{n=2} w_i f (x_i) 
     = f \left( - \frac{1}{\sqrt{3}} \right) + f \left( \frac{1}{\sqrt{3}} \right)
     = \frac{2}{3}
\end{align}
```
合っている．

$n = 3$とする．
```math
\begin{align}
    w_1
    &= \int_{-1}^{1} l_1 (x) \; dx
     = \frac{5}{9}
    \\
    w_2
    &= \int_{-1}^{1} l_2 (x) \; dx
     = \frac{8}{9}
    \\
    w_3
    &= \int_{-1}^{1} l_3 (x) \; dx
     = \frac{5}{9}
    \\
    I_{\text{GL}} [f; n=3]
    &= \sum_{i=1}^{n=3} w_i f (x_i) 
     = - \frac{5}{9} \frac{18}{5} + \frac{8}{9} 3
     = \frac{2}{3}
\end{align}
```
やはり合っている．

## 4次多項式
先の関数の次数を1だけ上げて，
```math
\begin{align}
    f (x)
    &= 7 x^4 - 8 x^3 - 3 x^2 + 3 x
\end{align}
```
を$[-1, 1]$で積分する．厳密解は$4 / 5 = 0.8$．

**台形則**
端点を集める．
```math
\begin{align}
    I_{\text{Trapz}} [f]
    &= 9 - 1
    \\
    &= 8
\end{align}
```
相変わらずpoor．

**Newton-Cotes**
$n = 3$として，$(x_1, x_2, x_3) = (-1, 0, 1)$を採用する．
```math
\begin{align}
    I_{\text{NC}} [f; n=3]
    &= \frac{1}{3} f \left( -1 \right)
        + \frac{4}{3} f \left( 0 \right)
        + \frac{1}{3} f \left( 1 \right)
     = \frac{3}{8} \approx 2.666
\end{align}
```
まだ足りない．

$n = 4$として，$(x_1, x_2, x_3, x_4) = (-1, -1/3, 1/3, 1)$を採用する．
```math
\begin{align}
    I_{\text{NC}} [f; n=4]
    &= \frac{1}{4} f \left( -1 \right)
        + \frac{3}{4} f \left( - \frac{1}{3} \right)
        + \frac{3}{4} f \left(   \frac{1}{3} \right)
        + \frac{1}{4} f \left( 1 \right)
     = \frac{44}{27} \approx 1.629 
\end{align}
```
まだ足りない．

$n = 5$として，$(x_1, x_2, x_3, x_4, x_5) = (-1, -1/2, 0, 1/2, 1)$を採用する．
```math
\begin{align}
    I_{\text{NC}} [f; n=5]
    &= \frac{7}{45} f \left( -1 \right)
        + \frac{32}{45} f \left( - \frac{1}{2} \right)
        + \frac{12}{45} f \left(   0           \right)
        + \frac{32}{45} f \left(   \frac{1}{2} \right)
        + \frac{7}{45} f \left( 1 \right)
     = \frac{4}{5}
\end{align}
```
足りた．

**Gauss-Legendre**
$n = 2$とする．
```math
\begin{align}
    I_{\text{GL}} [f; n=2]
    &= f \left( - \frac{1}{\sqrt{3}} \right) + f \left( \frac{1}{\sqrt{3}} \right)
     = - \frac{4}{9}
\end{align}
```
足りない．

$n = 3$とする．
```math
\begin{align}
    I_{\text{GL}} [f; n=3]
    &= \frac{5}{9} f \left( - \sqrt{\frac{3}{5}} \right)
        + \frac{8}{9} f \left( 0 \right)
        + \frac{5}{9} f \left( \sqrt{\frac{3}{5}} \right)
     = \frac{4}{5}
\end{align}
```
足りた．

今回は4次の項が落ちていないので，$n$点だけ与えられたとき，Newton-Cotesが$(n-1)$次以下の多項式を，Gauss-Legendreが$(2 n-1)$次以下の多項式を厳密に計算できる性質をrealizeできた．

## trigonometric
多項式でない
```math
\begin{align}
    I
    := \int_{0}^{\pi / 2} 
        \sin (x)
    \; dx
     = 1
\end{align}
```
を考える．Maclaurin展開を与えれば，（主観だが）5次程度の多項式で良く近似できる．

| $\text{func}$ | $\text{error}$ |
|:-:|:-:|
| ![fig (10).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/6d9e2e4a-547c-4f24-be71-a1e4c330e818.png) | ![err_trig.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/7312a6f4-b709-44d9-b217-53045e79e89d.png) |

[いまさら聞けない計算力学の定石](https://www.maruzen-publishing.co.jp/book/b10112229.html)で取り上げられている$I := \int_{0}^{\pi / 2} x \sin (x) \ dx = 1$に対しても，殆ど同じ結果が得られる．

## logarithmic
[いまさら聞けない計算力学の定石](https://www.maruzen-publishing.co.jp/book/b10112229.html)から，
```math
\begin{align}
    I
    := \int_{0}^{1}
        x \log (1 + x)
    \; dx
     = \frac{1}{4}
\end{align}
```
を考える．Maclaurin展開
```math
\begin{align}
    \log (1 + x)
    &= \sum_{n=1}^{\infty} \frac{(-1)^{n+1} x^{n}}{n}
    \\
    x \log (1 + x)
    &= \sum_{n=1}^{\infty} \frac{(-1)^{n+1} x^{n+1}}{n}
\end{align}
```
を使えば，

| $\text{func}$ | $\text{error}$ |
|:-:|:-:|
| ![fig (11).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/a3f62974-c754-4019-a60a-d019fa8599c5.png) | ![err_log.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/fa14b13b-6376-43a8-bd3a-941d1ab737c7.png) |

であり，相も変わらず強力である．

# 終わりに
Lagrange多項式が美しい．

| $n=2,\ p=1$ | $n=3,\ p=2$ | $n=4,\ p=3$ |
|:-:|:-:|:-:|
| ![lagrange_order1 (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/61943bc1-6e4c-4dde-bd56-bab9c0003ab4.png) | ![lagrange_order2 (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/f8e3ab61-b50e-4da6-8c9e-d5a3e9981924.png) | ![lagrange_order3 (1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/0ec4c851-afc3-4107-b044-c8ec8b3017cd.png) |

また，Legendre多項式も美しい．

| $P^{(0)}$ | $P^{(1)}$ | $P^{(2)}$ |
|:-:|:-:|:-:|
| ![legendre_deriv0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/2b5c53fe-9974-4392-8597-afa7346916b0.png) | ![legendre_deriv1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/7421fbfd-4e5f-4a30-88d7-20d5d8ed54c0.png) | ![legendre_deriv2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/35469580-36f1-4e86-bfe3-de7a6d44e85b.png) |