【勉強備忘録】Fortranで数値計算：ガウス求積をするmoduleを作ろう

# 概要

　Fortran90を用いて，ガウス求積(Gaussian quadrature)を行うmoduleを作成したので，その考え方や作成手順，スクリプトなどを記録するために本記事を作成しました．筆者はFortran自体学習中の言語であり，プログラムの速度や効率性はあまり考えられていないことを申し添えておきます．

### 更新履歴
2023/06/08 本記事を限定公開記事として投稿しました．
2023/09/09 内容を編集しました．
2023/11/26 本記事を公開設定にしました．

### 目次
1. [開発環境](#開発環境)
1. [ガウス求積の数式表現](#ガウス求積の数式表現)
1. [補足：Legendre多項式の1階微分の漸化式ver.の導出](#補足：Legendre多項式の1階微分の漸化式ver.の導出)
1. [moduleの設計](#moduleの設計)
1. [動作テスト](#動作テスト)
1. [参考文献](#参考文献)


# 開発環境
- OS: Windows11
- 言語: Fortran90/95

# ガウス求積の数式表現
次のように関数や文字を表します．
- $P_n(x)$：n次のLegendre多項式（漸化式） 
- $P_n^{'} (x)$：n次のLegendre多項式の1階微分（漸化式）※導出：
- $x_i(1\leq i \leq n)$：Gauss点（$P_n(x_i)=0$を満たす）
- $w_i(1\leq i \leq n)$：重み

このとき，上記の関数や文字は以下の式で定義されます．
```math
\begin{align}
P_{n}(x) &= \frac{1}{n}\left( (2n-1)xP_{n-1}(x)-(n-1)P_{n-2}(x) \right),
P_{0}(x) = 1, P_{1}(x) = x \\
P_{n}^{'}(x) &= \frac{1}{n}\left( (2n-1)xP_{n-1}^{'}(x)-(n-1)P_{n-2}^{'}(x) \right) + \frac{2n-1}{2}P_{n-1}(x) \\
w_i &= \frac{2}{ (1-{x_i}^2) { P_{n}^{'}(x) }^2 }\\
\int_a^b f(x) dx &= \frac{b-a}{2}\left( \sum_{i = 1}^n w_if\left( \frac{b-a}{2}x_i + \frac{a+b}{2} \right) \right)
\end{align}
```

### 補足：Legendre多項式の1階微分の漸化式ver.の導出
　n次のLegendre多項式の漸化式の両辺を$x$で微分して整理することで，n次のLegendre多項式の1階微分を漸化式の形で定義できあます．以下に簡易的に導出の流れを示しました．

まずは漸化式の両辺を微分します．
```math
\begin{align}
\frac{d}{dx}P_{n}(x) &= 
\frac{d}{dx}\left( \frac{1}{n}\left( (2n-1)xP_{n-1}(x)-(n-1)P_{n-2}(x) \right) \right)\\
P_{n}^{'}(x) &= \frac{1}{n}\left( (2n-1)\frac{d}{dx}(xP_{n-1}(x))-(n-1)\frac{d}{dx}P_{n-2}(x) \right)\\
\end{align}
```
積の微分に注意して微分と式変形を行います．

```math
\begin{align}
P_{n}^{'}(x) &= \frac{1}{n}\left( (2n-1)(P_{n-1}(x)+xP_{n-1}^{'}(x))-(n-1)P_{n-2}^{'}(x) \right)\\
P_{n}^{'}(x) &= \frac{1}{n}\left( (2n-1)xP_{n-1}^{'}(x)-(n-1)P_{n-2}^{'}(x) \right) + \frac{2n-1}{2}P_{n-1}(x)\\
\end{align}
```
ということで漸化式を得ました．ここで$n=0,n=1$のときは，

```math
\begin{align}
P_{0}^{'}(x) &= \frac{d}{dx}(1)=0\\
P_{1}^{'}(x) &= \frac{d}{dx}(x)=1\\
\end{align}
```
となるので，これで任意の$n$に対して再帰的な計算によって$P_{n}^{'}(x)$が求められます．

# moduleの設計
　今回のプログラムでは次のようなsubroutineを作成しました．

- legendre_polynomial：n次のルジャンドル多項式に，ある値$x$を代入したときの値を求める．（再帰可能）
- gauss_points：n次のルジャンドル多項式のGauss点を求める
- gauss_weight：n次のルジャンドル多項式の重みを求める
- diff_legendre_polynomial：n次のルジャンドル多項式の1階微分に，ある値$x$を代入したときの値を求める（再帰可能）

以下で各subroutineの仕様について詳しく説明します．

## legendre_polynomialの設計
　legendre_polynomial moduleはn次のルジャンドル多項式にある値$x$を代入したときの値を求めるものです．このmoduleではLegendre多項式の漸化式に則って再帰的に計算を繰り返すことで目的の値を導出するというアプローチをとっています．今回は採用しなかった別のアプローチとして，Legendre多項式の一般項を用いて，（再帰的ではない）計算を行うというものもあります．この場合，二項係数の計算が必要になりますので同一subroutine内または別のfunction/ subrouitneとして二項係数を計算するプログラムを定義した上で計算を行います．
　本moduleでは引数として次の3つの変数を渡しています．
- n：ルジャンドル多項式の次数
- x：ルジャンドル多項式に代入する値
- P：ルジャンドル多項式の計算結果
また内部の変数として次の2つの変数を定義しています．
- P_min1：n-1次のルジャンドル多項式の計算結果
- P_min2：n-2次のルジャンドル多項式の計算結果

```Fortran
recursive subroutine legendre_polynomial(n,x,P)
    integer(int32), intent(in) :: n
    real(real64), intent(inout) :: P,x
    real(real64) :: P_min1,P_min2

    !inductive definision
    !P_0(x) = 1, P_1(x) = x
    !P_{n}(x) = ((2n-1) x P_{n-1}(x) - (n-1) P_{n-2}(x)) / n
    
    if (n < 0) then
        print *, "error: n is too small"
    else if(n == 0) then
        P = 1.0d0
    else if(n == 1) then
        P = x
    else
        call legendre_polynomial(n-2, x, P_min2)
        call legendre_polynomial(n-1, x, P_min1)
        P = ((2.0d0*real(n)-1.0d0)* x * P_min1 - real(n-1) * P_min2)/real(n)
    endif
    
end subroutine legendre_polynomial
```
## gauss_pointsの設計
　全てのガウス点の座標を求めるsubroutineです．
ガウス点の判定方法は以下の通りです．まず領域を等分割し，その点におけるルジャンドル多項式のの値を求めます．次に隣り合う点同士の値の符号が異なる場合にはその間にガウス点が存在していることが言えるので，中点をガウス点として代入します．この手法のデメリットとしてルジャンドル多項式そのものの値が0になってしまう場合にその点がガウス点として扱われずに代入されないという事象が発生してしまうという点があります．本プログラムではそれを回避する直接的な手法は導入しませんでしたが，計算時にガウス点そのものを出力してユーザー側が確認をするという方法をとりました．この点を改善することは今後の課題です．
　またガウス点の性質としてルジャンドル多項式の次数nが奇数であればその点の中に0をかならず含むという性質を用いて，他の点の計算終了後にnが奇数の場合には配列gaussに0を代入しています．
- n：ルジャンドル多項式の次数，すなわちガウス点の数．
- gauss：ガウス点（配列）
- P：ルジャンドル多項式の値
- nsize：分割の細かさ
- i, j：ループ用の変数

```Fortran
subroutine find_gauss_points(n,gauss)
    integer(int32), intent(in):: n
    real(real64), intent(inout) :: gauss(:)
    real(real64), allocatable :: P(:,:) !=(x, P(x))
    integer(int32) :: nsize,i,j
    ! n: number of Gauss points
    ! gauss(:): Gauss points

    nsize = 5000 !number of discrete points
    allocate(P(nsize,2))
    P = 0.0d0
    j = 1
    
    P(1,1) = -1.0d0
    call legendre_poly(n, P(1,1), P(1,2))

    do i = 2, nsize
        P(i,1) =  -1.0d0 + 2.0d0/dble(nsize-1) * dble(i-1)
        call legendre_poly(n, P(i,1), P(i,2))

        if(P(i-1,2)*P(i,2) < 0) then
            gauss(j) = (P(i-1,1) + P(i,1))/2.0d0
            j = j + 1
        endif
    enddo
    !when n is odd number, 0 must be the gauss point
    if(mod(n,2) == 1) then
        gauss((n+1)/2) = 0.0d0
    endif

end subroutine find_gauss_points
```

## gauss_weightの設計
　このsubroutineでは各ガウス点における重みwを計算するプログラムです．式に直接値を代入することでその重みを求めています．ただしここで，n次のルジャンドル多項式の1階微分にガウス点の値を代入した結果が欲しいので，diff_legendre_polynomialというsubroutineを用意して計算を別に実行しています．
- n：ルジャンドル多項式の次数，すなわちガウス点の数．
- w(:)：重み
- gauss(:)：ガウス点
- i：ループ用の変数
- Pd：n次のLegendre多項式の1階微分
```Fortran
subroutine gauss_weight(n,gauss,w)
        integer(int32), intent(in) :: n
        real(real64), intent(inout) :: w(:)
        real(real64), intent(inout) :: gauss(:)
        integer(int32) :: i
        real(real64) :: Pd

        do i = 1, n
            !weights expressed by (n-1)order legendre polynomial
            call diff_legendre_poly(n, gauss(i), Pd)
            w(i) = 2.0d0 / ((1.0d0-gauss(i)*gauss(i)) * (Pd*Pd))
        enddo

    end subroutine gauss_weight
```

## diff_legendre_polynomialの設計
 このプログラムでは[上記](###補足：Legendre多項式の1階微分の漸化式ver.の導出)で示した，n次のルジャンドル多項式の1階微分を求める漸化式を用いて，その具体的な値の計算を行います．設計としては単に数式をそのままプログラムで表現しています．

- n：ルジャンドル多項式の次数，すなわちガウス点の数．
- x
- Pd：n次のLegendre多項式の1階微分
- P_nmin1：n-1次のLegendre多項式
- Pd_nmin1：n-1次のLegendre多項式の1階微分
- Pd_nmin2：n-2次のLegendre多項式の1階微分


```Fortran
recursive subroutine diff_legendre_poly(n,x,Pd)
    integer(int32), intent(in) :: n
    real(real64), intent(inout) :: x
    real(real64) :: Pd, P_nmin1, Pd_nmin1, Pd_nmin2
    
    if(n == 0) then
        Pd = 0
    else if(n == 1) then
        Pd = 1
    else
        call legendre_poly(n-1,x,P_nmin1)
        call diff_legendre_poly(n-1,x,Pd_nmin1)
        call diff_legendre_poly(n-2,x,Pd_nmin2)
        Pd = 1.0d0/real(n) * (real(2*n-1) * x * Pd_nmin1 &
                            - real(n-1)       * Pd_nmin2)&
                            + real(2*n-1)/real(n) * P_nmin1
    endif

    end subroutine diff_legendre_poly
```

# 動作テスト
　では以上のsubroutineを用いて動作確認をしてみます．なお，ガウス点を求める際，領域[-1,1]の分割数は一律で5,000としました．

　それでは次のような２つの積分を実行してみます．
```math
    \int_1^{10} e^xdx = e^{10}-e^1 \simeq 22023.74751 \cdots \textrm{①}
```
```math
    \int_1^{10} \ln{x} =10\ln{(10)} -9 \simeq 14.02585092 \cdots \textrm{②}
```
なお右辺の値は厳密解を適当な位置でで四捨五入した値で，本テストではこの値との相対誤差で性能を検証します．またルジャンドル多項式の次数は1~11とします．

|ルジャンドル多項式の次数|積分①の結果|①の相対誤差	|積分②の結果	|②の相対誤差|
|---|---|---|---|---|
|1	|2202.22739     |0.900006691	|15.34273283	|0.093889627|
|2	|14877.14695	|0.324495209	|14.20751432	|0.012952042|
|3	|20966.77799	|0.047992265	|14.06012466	|0.002443612|
|4	|21944.42311	|0.003601767	|14.03704851	|0.000798354|
|5	|22019.72513	|0.000182638	|14.0267904	    |6.6982E-05|
|6	|22010.43359	|0.000604526	|14.02171324	|0.000295004|
|7	|22038.1578	    |0.000654307	|14.03313185	|0.000519108|
|8	|22041.61227	|0.000811159	|14.02944299	|0.000256103|
|9	|22014.88973	|0.000402192	|14.02416852	|0.00011995|
|10	|22045.9734	    |0.001009179	|14.03160411	|0.000410185|
|11	|22031.80586	|0.000365894	|14.02884456	|0.000213437|

まず①の積分について，その値と誤差の変化を見てみます．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2261537/a079725a-5b67-cc38-9984-ed707d37b176.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2261537/82b62ded-f345-6a87-7e24-8c91648861ba.png)
n=5以降で比較的変化が緩やかになっていることがわかります．ただ，n=5以降の値の変化に注目すると面白い現象が見えました．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2261537/417f5f0e-8d09-cdc5-f2e4-edc0c9e094a3.png)
グラフから分かるように，nの増加に伴って誤差が単調に減少するわけではなく，増減が見られます．また，nが奇数の時に比較的誤差の小さい結果が得られていることもわかります．これは，ガウス点を求めるsubroutineの仕様である，nが奇数のときには強制的に0をガウス点の座標として追加するという効果が表れているものだと考えられます．


同様に②の積分結果についても見てみます．こちらについてはn=3以降で変化が緩やかになっています．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2261537/4040c720-3597-9ddf-b062-b941406c98f2.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2261537/8b8bd946-a9e9-4033-03b7-08f8f11bd132.png)
こちらについても①の積分と同様にn=3以降の変化に注目すると，相対誤差は単調に減少せず奇数の場合に精度が高くなる場合が存在していることが分かりました．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2261537/1f6b5f9b-f1f1-83bb-3053-f84b63eae8c0.png)




# 参考文献
　ガウス求積の数学的な理論に関しては以下のWikipediaのページを参考にしました．（日本語版と英語版で記載内容が若干異なるので両方掲載します）
- [ガウス求積 - Wikipedia](https://ja.wikipedia.org/wiki/%E3%82%AC%E3%82%A6%E3%82%B9%E6%B1%82%E7%A9%8D)
- [Gaussian quadrature - Wikipedia](https://en.wikipedia.org/wiki/Gaussian_quadrature)