# 積分方法第一部
{%hackmd @themes/orangeheart %}
<style>
.likecoin-button {
position: relative;
width: 100%;
max-width: 485px;
max-height: 240px;
margin: 0 auto;
}
.likecoin-button > div {
padding-top: 49.48454%;
}
.likecoin-button > iframe {
position: absolute;
top: 0;
left: 0;
width: 100%;
height: 100%;
}
</style>
###### tags: `Numerical Methods`
在統計學或經濟學中，我們常常會碰到以下的積分問題：
$$
\begin{aligned}
I = \int_a^b p(x) f(x) dx,
\end{aligned}
$$
其中 $p(x)$ 是一個非負的**加權函數(weight function)**。例如，根據不自覺統計學家法則，假設給定一個機率密度函數$p(x)$，若要計算他的期望值，我們便可以利用上述的方法計算出 $I = \mathbb{E}[f(x)]$。如果說這個加權函數等於$1$，那麼積分的意義就變成是計算從下界$a$到上界$b$與該曲線所圍成的面積。
## 定積分的近似值
給定一個區間 $[a,b] = [2, 20]$ 以及一個函數$g(x) = \frac{1}{1+x}$，而我們想要計算在給定的區間時曲線下的面積，即
$$
\begin{aligned}
I = \int_a^b g(x) d x = \int_2^{20} \frac{1}{1+x} dx.
\end{aligned}
$$
而根據計算機所算出的閉合解(close form)，得到 $I = \log(7) \approx 1.9459\cdots$。我們可以在 $\texttt{Julia}$ 畫個圖。
根據我們在微積分課程學到的方法（大概率都是教黎曼和），我們可以在給定的區間內，將曲線下的面積切成好多個矩形並加總，當切成越多塊其就會越近似我們的答案。回到上面的例子，我們可以把區間分成10塊，即 $[2,4,6,8,10,12,14,16,18,20]$，接著可以在$[2,4]$、$[4,6]$$\cdots$ $[18,20]$找出節點(node)，就可以進行加總運算。
### 矩形加總法(Forward Rectangular Rule)
我們利用最左邊的端點當作區間，找出相對應的矩形進行計算。第一個矩形為$(4-2)g(2) = (4-2)\frac{1}{1+2}$，第二個則是$(6-4)g(4) = (6-4)\frac{1}{1+4}$ $\cdots$。因此我們可以得到
$$
\begin{aligned}
I = \int_2^{20}\frac{1}{1+x}dx \approx (4-2)\frac{1}{1+2} + (6-4)\frac{1}{1+4} + \ldots,
\end{aligned}
$$
我們可以推廣上述的計算方式：
$$
\begin{aligned}
I = \int_a^{b} g(x) dx \approx & \sum_{i=1}^n \omega_i g(\xi_i) \\
= & \frac{b-a}{n} \sum_{i=1}^n g(\xi_i)
\end{aligned}
$$
在我們上述的例子中，區間$[a, b] = [2, 20]$，$n = 9$，$\omega_i = \frac{b-a}{n} = 2$則是一個矩形（區間）的長度，而 $\{\xi_i\} = \{2,4,6,8,10,12,14,16,18 \}$。
```jula
using Interact, WebIO, StatsPlots
legend1 = LaTeXStrings.LaTeXString("\$\\frac{1}{1+x}\$")
g(x) = 1/(1+x)
true_value = 1.9459101490553133051053527434432
w_i = 3:2:19
g_i = @.g(w_i-1)
bar(w_i, g_i, xticks=0:2:22, bar_width=2, fillalpha = 0.35, label="")
plot!(g,xlim=[0,22], ylim=[0.0,0.4], xlabel="x", ylabel="g(x)", title="Forward Rectangular Rule",
label=legend1, size=(450,280)) |> display
@manipulate for n in (2, 3, 6, 9, 18, 36, 90, 180)
# @manipulate for n in 1:90
step = (20-2)/n;
w_i = 2+step/2:step:20-step/2;
g_i = g.(w_i .- step/2);
res1 = step*sum(g_i)
err1 = res1 - true_value
weights = fill(step, n);
nodes = collect(w_i .- step/2);
rule = collect(zip(weights, nodes));
@show weights
@show nodes
@show rule
bar(w_i, g_i, xticks=0:2:22, bar_width=step, fillalpha = 0.35, label="area=$(res1)\n err=$(err1)")
plot!(g, xlim=[0, 22], ylim=[0.0, 0.4], xlabel="x", ylabel="g(x)", title="Forward Rectangular Rule",
label="", size=(650,380))
end
```
### 中點法(Midpoint Rule
從上圖中我們可以看到，矩形與曲線相切時會出現誤差，而越接近最初始的矩形則誤差越大。而中點法可以試圖解決這個問題：我們找到一個區間的兩個點，並找出其中點。承接上面的例子，第一個矩形的計算方式改為$(4-2)*g((4+2)/2) = (4-2)*\frac{1}{1+(4+2)/2}$，其中$\frac{4+2}{2}$ 代表 $2$ 與 $4$ 的中點。同樣地，我們可以推廣它
$$
\begin{aligned}
I = \int_a^b g(x) dx \approx & \sum_{i=1}^n \omega_i g(\tilde{\xi}_i) \\
= & \frac{b-a}{n} \sum_{i=1}^n g(\tilde{\xi}_i),
\end{aligned}
$$
在我們的例子中，$\{\tilde{\xi}_i\} =\{(\xi_i + \xi_{i+1})/2\} = \{3,5,7,9,11,13,15,17,19 \}$
```julia
g2(x) = 1/(1+x)
true_value = 1.9459101490553133051053527434432
x_i = 3:2:19
g_i = g2.(x_i)
bar(x_i, g_i, xticks=1:2:21, bar_width=2, fillalpha = 0.35, label="")
plot!(g2, xlim=[0,22], ylim=[0.0,0.4], title="Midpoint Rule", size=(450, 280),
label=legend1, xlabel="x", ylabel="g2(x)") |> display
@manipulate for n in (2, 3, 6, 9, 18, 36, 90, 180)
step = (20-2)/n
x_i = 2+step/2:step:20-step/2
g_i = g2.(x_i)
res2 = step * sum(g_i)
err2 = res2 - true_value
weights = fill(step, n);
nodes = collect(x_i);
rule = collect(zip(weights, nodes));
@show weights
@show nodes
@show rule
bar(x_i, g_i, xticks=1:1:21, bar_width=step, fillalpha = 0.35, label="area=$(res2)\n err=$(err2)")
plot!(g2, xlim=[0,22], ylim=[0.0,0.4], title="Midpoint Rule", size=(650, 380),
label="", xlabel="x", ylabel="g2(x)")
end
```
### 梯形法(Trapezoidal Rule)
矩形加總法與中點法都是利用單點找出矩形，然後進行加總。當然，我們也可以找出兩點並加總，而梯形法正是其中一個例子：找出端點與區間並繪出梯形、加總。
$$
\begin{aligned}
I = \int_a^b g(x) dx \approx & \sum_{i=1}^n \omega_i \frac{ g(\xi_i) + g(\xi_{i+1})}{2} \\
= & \frac{b-a}{n} \sum_{i=1}^n \frac{ g(\xi_i) + g(\xi_{i+1})}{2},
\end{aligned}
$$
```julia
g3(x) = 1/(1+x)
true_value = 1.9459101490553133051053527434432
w_i = 2:2:20
g_i = @.g3(w_i)
plot(w_i, zeros(length(w_i)), fillrange=g_i, fillalpha = 0.35, label="")
plot!(w_i, g_i, xticks=0:2:22, shape=:circle, markersize=2, label="", seriestype = :scatter)
plot!(g3, xlim=[0, 22], ylim=[0.0, 0.4], title="Trapezoidal Rule", size=(450, 280),
label=legend1, xlabel="x", ylabel="g(x)") |> display
@manipulate for n in (2, 3, 6, 9, 18, 36, 90, 180)
step = (20-2)/n
w_i = 2:step:20
g_i = @.g3(w_i)
g_i2 = [g_i[j] for j in 2:n+1]
res3 = (step/2)*sum(g_i[1:n] .+ g_i2)
err3 = res3 - true_value
plot(w_i, zeros(length(w_i)), fillrange=g_i, fillalpha = 0.35, label="area=$(res3)\n err=$(err3) ")
plot!(w_i, g_i, xticks=0:2:22, shape=:circle, markersize=2, label="", seriestype = :scatter)
plot!(g3, xlim=[0, 22], ylim=[0.0, 0.4], title="Trapezoidal Rule",
label="", xlabel="x", ylabel="g(x)")
end
```
:::warning
經過計算，我們可以發現矩形加總計算出來的積分數值，其誤差為梯形法的一半！
:::
既然我們觀察到上面的現象，或許我們把這些方法**加總起來**就可以弭平誤差了。事實上，這個想法便是[辛普森法](https://en.wikipedia.org/wiki/Simpson%27s_rule)背後的邏輯。
### 辛普森法(Simpson's Rule)
$$
\begin{aligned}
I = \int_a^b f(x) dx \approx \frac{b-a}{6n} \sum_{i=1}^n \left[ f(\xi_i) + 4f\left(\frac{\xi_i + \xi_{i+1}}{2}\right) + f(\xi_{i+1}) \right].
\end{aligned}
$$
它利用了三個點： $\xi_i$、$(\xi_i + \xi_{i+1})/2$ 與 $\xi_{i+1}$，權重則為 $\frac{b-a}{6n}$、$\frac{2(b-a)}{3n}$ 與 $\frac{b-a}{6n}$。有趣的是，辛普森法的採用的 $n$ 個節點 $I^{S(n)}$ 是梯形法($T$)與中點法($M$)的**加權平均**（因兩者均使用 $\frac{n}{2}$ 個節點），即：
$$
\begin{aligned}
I^{S(n)} = \frac{1}{3}I^{T\left(\frac{n}{2}\right)} + \frac{2}{3}I^{M\left(\frac{n}{2}\right)}.
\end{aligned}
$$
```julia
using Interact, WebIO, StatsPlots
# Simpson's Rule
g4(x) = 1/(1+x)
true_value = 1.9459101490553133051053527434432
@manipulate for n in (2, 3, 6, 9, 18, 36, 90, 180) # number of segments
a, b= 2, 20
step = (b-a)/n # length of one segment
x = [a:step:b;] # points that divide the segments
width = [(b-a)/(6n), 2(b-a)/(3n), (b-a)/(6n)] # coefficients of the 3 points = width of bars
ω = repeat(width, outer=n) # width of bars, spread to all segments
start = repeat(x[1:n], inner=3) # start of each segment (each seg has 3 bars); don't need the last elem of x
bar_dist = repeat([0.5*width[1], width[1]+0.5*width[2], width[1]+width[2]+0.5*width[3]], outer=n) # mid of bar measured from the beginning of segment
bar_mid = start .+ bar_dist # mid of each bar, for plotting purposes
ξ_dist = repeat([0, width[1]+0.5*width[2], width[1]+width[2]+width[3]], outer=n) # ξ measured from the beginning of segment
ξ = start .+ ξ_dist # the ξ in the formula
g_height = g4.(ξ) # the height of the bar
res = sum(ω .* g4.(ξ)) # the area = solution
err = res - true_value
rule = collect(zip(ω, ξ));
@show ω
@show ξ
@show rule
bar(bar_mid, g_height, xticks=1:1:21, bar_width=ω, fillalpha=0.35, label="area = $res\n err = $err") # bar plot
plot!(g4, xlim=[0, 22], ylim=[0.0, 0.4], title="Simpson's Rule", label="", xlabel="x", ylabel="g(x)") # plot g4
scatter!(ξ, g4.(ξ), label=false) # highlight the dividing points
end
```
## 數值積分法(Quadrature Methods)與多項式積分法(Polynomial Integration Methods)
上述幾個直觀的方法被歸類為數值積分法：利用**權重**與**節點**進行函數積分的估計。當然，我們也可以利用其他的方法，目的正是為了將誤差變小，或是對於其他不同類型的函數進行積分近似。我們可以將上述的幾個方法進行小結，數值積分法的想法是將函數與權重函數進行相乘積，並透過權重與節點進行近似估計：
$$
\begin{aligned}
I = \int_a^b p(x)f(x) dx \approx \sum_{i=1}^n \omega_i f(\xi_i),
\end{aligned}
$$
不過，必須注意到以下幾個要點：
- $p(x)$ 與 $\omega_i$ 是兩個完全不同的東西，雖然都叫做權重，不過前者指涉的通常是一個正函數，而後者則是根據函數計算出的權重。
- $p(x)$ 事實上是被積函數(integrand)的一部分，不過不會出現在等式的右邊(RHS)。
- $\{\omega_i, \xi_i\}$， $i=1,\cdots,n$ 的這個性質與被積函數(integrand)**相互獨立**。
我們在此可以得到一個結論：數值積分法正是一個不斷尋找 $\{\omega_i, \xi_i\}$ 的方法，而不同取值的方式會對於其近似程度（誤差）與效率有很大的關係[^1]。針對不同的 $\{\omega_i, \xi_i\}$ 有各種取值的方法，而這些方法則是依照被積函數而有不同。而在實務上最常被使用的方法是**高斯積分法(Gauss Quadrature Methods)**。在上述的例子中，我們的上下界為 $[a, b ] = [2, 20]$，$p(x) = 1$，$f(x) = \frac{1}{1+x}$，$\omega_i = 2$，而$\xi_i$則是端點(endpoints)或是中點(midpoint)。
### 高斯積分法
高斯積分法基本上是處理不同類型的積分，並對於給定不同的上下界有專門的處理方式。我們可以給定**節點**數，產生$\{\omega_i, \xi_i\}$， $i=1,\cdots,n$。而常見的高斯積分法有以下四種類型，可以透過下表整理呈現[^2]：
| 積分法 | 值域 | $$p(x)$$ | 近似值 |
| :- | :- | :-: | :-: |
|Gauss-Legendre |$$[-1,1]$$ |$$1$$ |$$\int_{-1}^1 f(x) \approx \sum_{i=1}^n \omega_i f(\xi_i)$$ |
|Gauss-Chebyshev |$$[-1,1]$$| $$(1-x^2)^{-\frac{1}{2}}$$ |$$\int_{-1}^1 p(x)f(x) \approx \sum_{i=1}^n \omega_i f(\xi_i)$$ |
|Gauss-Hermite |$$(-\infty,\infty)$$| $$e^{-x^2} $$| $$\int_{-\infty}^\infty p(x)f(x) \approx \sum_{i=1}^n \omega_i f(\xi_i)$$ |
|Gauss-Laguerre |$$[0,\infty)$$ |$$e^{-x}$$ |$$\int_0^\infty p(x)f(x) \approx \sum_{i=1}^n \omega_i f(\xi_i)$$ |
那麼，我們應該用何種積分法進行積分運算呢？事實上，根據上表中的值域我們不難看出，針對不同的上下界，我們就要使用不同的積分方法。我們以 Gauss-Legendre 積分法為例：
- 給定節點數，根據 Legendre 多項式：
$$
(n+1)P_{n+1}(\xi) = (2n+1)\xi P_n(\xi) - n P_{n-1}(\xi)
$$
透過遞迴(recursively)的方式，可得到 $\{\xi_i\}$。其中 $P_0(\xi)=1$， $P_1(\xi)=\xi$。
- 而權重 $\{\omega_i\}$ 的計算方式可以藉由下列式子得到：
$$
\omega_{i}={\frac {2}{\left(1-\xi_{i}^{2}\right)\left[P'_{n}(\xi_{i})\right]^{2}}}
$$
在實務上有各種不同計算節點數與權重的方式，例如取得一個特定矩陣的特徵值(eigenvalues)與特徵向量(eigenvectors)。[^3]
### 積分變數變換
令 $x= \rho(t)$，則
$$
\begin{aligned}
\int_a^b g(x) dx = \int_{\rho^{-1}(a)}^{\rho^{-1}(b)} g(\rho(t)) \rho'(t) dt,
\end{aligned}
$$
其中的 $\rho^{-1}(t)$ 我們稱之為 Jacobian 項。我們簡單的來說明一下何謂 Jacobian項，以及其與梯度(gradient)及Hessian矩陣的差別。
- **梯度(Gradient)**：梯度本質上是方向導數，其是一個由函數 $f$ 從 $\mathbb{R}^n$ 映射到[^4] $\mathbb{R}$ 組成的向量。例如給定以下的 $f$：
$$f(x,y) = x^2 + xy$$
其梯度為
$$(2x+y, x)$$
- **Hessian 矩陣**：其為一個函數的二階偏導數構成的矩陣。同樣地給定$f$：
$$f(x,y) = x^2 +xy$$
則 Hessian 矩陣為
\begin{bmatrix}
2, 1 \\
1, 0
\end{bmatrix}
- **Jacobian 項**： 雅可比矩阵是函数的一阶偏导数以一定方式排列成的矩阵 Jacobian 項是由一個函數的一階偏導數以一定方式排列的矩陣。如
$$f(x,y) = (x^2 + xy, x^3 - 2xy)$$
其 Jacobian 項為
\begin{bmatrix}
2x+y, x\\
3x^2-2y, -2x
\end{bmatrix}
以下給出有用的 Jacobian 項與轉換法[^2]：
x 值域 | 轉換規則 | t 值域 | Jacobian 項|
:---: | :---: | :--- | ---
$$[a, b]$$ | $$x = \frac{a+b}{2} + \frac{b-a}{2}t$$ | $$[-1,1]$$ | $$\frac{b-a}{2}$$
$$[-\infty, \infty]$$ | $$x = \frac{t}{1-t^2}$$ | $$[-1,1]$$ | $$\frac{t^2 + 1}{(t^2 -1)^2}$$
$$[a, \infty]$$ | $$x = a + \frac{1+t}{1-t}$$ | $$[-1,1]$$ | $$\frac{2}{(t-1)^2}$$
$$[-\infty, b]$$ | $$x = b - \frac{1-t}{1+t}$$ | $$[-1,1]$$ | $$\frac{2}{(t+1)^2}$$
### Gauss-Legendre
Gauss-Legendre 積分法是最常見的數值積分方法，通常被積函數的形式都會是：
$$
\begin{aligned}
I = \int_a^b g(x) d x.
\end{aligned}
$$
因此我們要將其上下界轉換為 $[-1, 1]$。回到上面的例子，我們如果要計算
$$
\begin{aligned}
I = \int_a^b g(x) dx = \int^{20}_{2} \frac{1}{1+x} dx
\end{aligned}
$$
首先，其上下界為 $[2, 20]$，參照上方的轉換規則，可以得到
$$
\begin{aligned}
& x = \rho(t) = \frac{a+b}{2} + \frac{b-a}{2}t
\end{aligned}
$$
要檢驗的方式就是將 $t$ 的上下界代入，如果 $t$ 代 $-1$，就會得到 $a$；如果 $t$ 代入 $1$，就會得到 $b$。接著我們要求出 Jacobian 項，作法就是對轉換規則取一次導函數，我們會得到
$$
\rho'(t) =\frac{\partial \rho(t)}{ \partial t} = \frac{b-a}{2}
$$
最後，只要將原本函數中的 $x$ 替換為轉換規則，再乘上 Jacobian 項即可。在 $\verb|Julia|$ 中，我們可以利用 `FastGaussQuadrature` 這個套件產生節點，接著再進行一些基本的設定就可以計算積分了。
$$
\begin{aligned}
I = \int_a^b g(x) dx & = \int_{-1}^1 g\left( \frac{1}{2}[ (b+a) + (b-a)t ] \right) \frac{b-a}{2}dt.
\end{aligned}
$$
```julia
## 引入套件
using FastGaussQuadrature
## 設定積分上下界
a = 2; b=20
## 設定節點數
n=30
## 設定被積函數
g(x) = 1/(1+x)
## 設定轉換規則與 Jacobian 項
trans_rule(t) = (a + b)/2 + t*(b-a)/2
jacobian = (b-a)/2
## 重新設定積分形式
f(t) = g( trans_rule(t) ) * jacobian
## 產生 x 與節點並計算
xi, wi = gausslegendre(n)
answer = sum(wi .* f.(xi))
## 1.9459101490553121
```
我們也可以用另一個套件 `QuadGK`進行檢驗，看看我們剛剛的方法是否正確。[^5]
```julia
using QuadGK
quadgk(x-> 1/(1+x), 2, 20, rtol=1e-8)
## 1.9459101490553135
```
可以看到我們利用 Gauss-Legendre 積分的結果與實際值相差不大。
[^1]: 從上述的矩形加總法到中點法就可看出，當取節點越少時，其誤差便會稍稍的增加。
[^2]:由國立臺灣大學經濟學系[王泓仁](https://homepage.ntu.edu.tw/~wangh/)教授整理。
[^3]: *Mathematics of Computation*, vol. 35, no. 151, 1980, https://doi.org/10.1090/mcom/1980-35-151.
[^4]: 數學上我們將**映射**到記作：$\mapsto$。
[^5]: `rtol` 是該套件中容忍相對誤差的指令。
<div class="likecoin-embed likecoin-button">
<div></div>
<iframe scrolling="no" frameborder="0" src="https://button.like.co/in/embed/xiaolong70701/button?referrer=hackmd.io"></iframe>
</div>
Sign in with Wallet
Connect another wallet