多次元の求積

# 初めに
[MC/QMC](https://qiita.com/ShotaDeguchi/items/4b2121385543e1c142fc)と[quad rules](https://qiita.com/ShotaDeguchi/items/2b78da4c500fb5f39301)の続き．何度味わっても飽きない．

Gauss-Legendreは多次元で折り畳んで重なるが，速度は如何ほどか？

# Methods
## Quadrature
次を計算したい．
```math
\begin{align}
    I \left[ f (\boldsymbol{x}), D \right]
    := \int_{\boldsymbol{x} \in D}
        f (\boldsymbol{x})
    \; d\boldsymbol{x}
    &\approx \sum_{\boldsymbol{x}_i \in D}
        w_i
        f (\boldsymbol{x}_i)
\end{align}
```

## Monte-Carlo
次をMonte-Carlo求積と呼ぶ．
```math
\begin{align}
    I_{\text{MC}}
    &= \frac{|D|}{N} \sum_{i=1}^{N}
        f (\boldsymbol{x}_i)
    \\
    \mathbb{E} \left[ I_{\text{MC}} \right]
    &= \frac{|D|}{N} \sum_{i=1}^{N}
        \mathbb{E}_{\boldsymbol{x}_i \sim \mathcal{U} (D)} \left[
            f (\boldsymbol{x}_i)
        \right]
    \\
    &= \int_{\boldsymbol{x} \in \Omega}
        f (\boldsymbol{x})
    \; d\boldsymbol{x}
    = I
\end{align}
```
これは等重み求積であり，明らかにpartition of unityを満たす．$\boldsymbol{x}_i$には，所謂Monte Carloから選ばれるものと，quasi-Monte Carloから選ばれるものを採用する．

## Gauss-Legendre
次をGauss-Legendre求積と呼ぶ．
```math
\begin{align}
    I_{\text{GL}}
    &= \sum_{i=1}^{N}
        w_i
        f (\boldsymbol{x}_i)
    \\
    w_i
    &= \int_{-1}^{1} l_i (\boldsymbol{x}) \; d\boldsymbol{x}
    \\
    l_i (\boldsymbol{x})
    &= \prod_{j = 1; j \neq i}^{n} \frac{\boldsymbol{x} - \boldsymbol{x}_j}{\boldsymbol{x}_i - \boldsymbol{x}_j}
    \\
    P_n (\boldsymbol{x}_i)
    &= 0
\end{align}
```
これは重み付き求積であり，やはりpartition of unityを満たす．なお，$l_i$はLagrange多項式の基底，$P_n$はLegendre多項式である．

# Results
## 2D
### Test problem 1-2D
以下を考える．
```math
\begin{align}
    \int_{[0, \pi/2]^2}
        \sin (x_1) 
        \sin (x_2) 
    \; d\boldsymbol{x}
    &= 1
\end{align}
```

| $f$ | $\epsilon$ |
|:-:|:-:|
| ![fig_trig_2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/e80220fb-0192-48e0-932a-94a8f0170288.png) | ![err_trig_2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/908d09ea-a5ba-4fdc-b8a1-11955079edfa.png) |

Maclaurin展開により，

| $\text{poly of order}$ | $2$ | $4$ | $6$ | $8$ |
|:-:|:-:|:-:|:-:|:-:|
| - | ![fig_trig_2d_approx_p2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/75c93665-8bcc-42ef-a008-7c249375f109.png) | ![fig_trig_2d_approx_p4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/23d0232a-1fbc-4d6e-a37c-17e6286edb4c.png) | ![fig_trig_2d_approx_p6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/e7b40075-9499-4bc0-a34d-29eaaad54e4c.png) | ![fig_trig_2d_approx_p8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/3ae2b762-2d9a-4125-8e48-4a139941667f.png) |

などによって近似していることが分かる．

### Test problem 2-2D
以下を考える．
```math
\begin{align}
    \int_{[0, 1]^2}
        x_1 \log (x_1) \
        x_2 \log (x_2) 
    \; d\boldsymbol{x}
    &= \frac{1}{16}
\end{align}
```

| $f$ | $\epsilon$ |
|:-:|:-:|
| ![fig_log_2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/c37074a7-c5c5-4f91-9c99-3a01353d3af8.png) | ![err_log_2d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/d3b00d5e-997e-4b63-b51e-4d86375260b9.png) |

やはりMaclaurin展開により，

| $\text{poly of order}$ | $\text{4}$ | $\text{5}$ | $\text{6}$ | $\text{7}$ |
|:-:|:-:|:-:|:-:|:-:|
| - | ![fig_log_2d_approx_p4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/5a12a9de-8b3e-4ce0-a929-0996138c122b.png) | ![fig_log_2d_approx_p5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/49e2d403-0013-4385-8e02-ed724c7023c3.png) | ![fig_log_2d_approx_p6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/28d79601-a311-4dd1-96d0-beaf55f4efe5.png) | ![fig_log_2d_approx_p7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/3fce2ac6-04f4-4983-8368-0af47350ed65.png) |

を以て近似している．

## 3D
### Test problem 1-3D
```math
\begin{align}
    \int_{[0, \pi/2]^3}
        \sin (x_1) 
        \sin (x_2) 
        \sin (x_3) 
    \; d\boldsymbol{x}
    &= 1
\end{align}
```

| $f$ | $\epsilon$ |
|:-:|:-:|
| ![fig_trig_3d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/70e3dc6b-5938-4d33-879b-9ab81b34f1fc.png) | ![err_trig_3d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/cc58de80-5d95-4b77-b7e8-c215ac40733a.png) |

### Test problem 2-3D
```math
\begin{align}
    \int_{[0, 1]^3}
        x_1 \log (x_1) \
        x_2 \log (x_2) \
        x_3 \log (x_3) 
    \; d\boldsymbol{x}
    &= \frac{1}{64}
\end{align}
```

| $f$ | $\epsilon$ |
|:-:|:-:|
| ![fig_log_3d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/6114a89b-8b98-472a-82e0-665d6ad05a50.png) | ![err_log_3d.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/6bcbbf0e-10da-459d-ae76-a0f82c559b2f.png) |

次元が大きくなると速度が落ちる．直感に一致する結果である．しかし，Gauss-LegendreはMCやQMCなど眼中に無いほどに速い．

# 終わりに
正直に言うと，Gauss-LegendreはGridと同様に重なってゆくから，余り期待していなかった．勿論，強力であることは認めていたし，必要以上に失望していたわけではないのだが．．．予想を遥かに上回るほどに凄まじい結果であった，ということである．

多次元程度ではGauss-Legendreの独壇場のようである．高次元ではMonte Carloにスポットライトが当たるだろうか．

| $\text{ID}$ | $\text{Quadrature nodes}$ |
|:-:|:-:|
| $\text{Test problem 1}$ | ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/2cfde1b3-9fb1-4cdb-9671-0d3b325b3e7b.png) |
| $\text{Test problem 2}$ | ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/ee9ae6e4-d3e0-44da-a741-ec7cdacedeca.png) |