概要

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

更新履歴

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

開発環境

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)$：重み

このとき，上記の関数や文字は以下の式で定義されます．

\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階微分を漸化式の形で定義できあます．以下に簡易的に導出の流れを示しました．

まずは漸化式の両辺を微分します．

\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}

積の微分に注意して微分と式変形を行います．

\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$のときは，

\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次のルジャンドル多項式の計算結果

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：ループ用の変数

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階微分

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の設計

このプログラムでは上記で示した，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階微分

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としました．

　それでは次のような２つの積分を実行してみます．

    \int_1^{10} e^xdx = e^{10}-e^1 \simeq 22023.74751 \cdots \textrm{①}

    \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

まず①の積分について，その値と誤差の変化を見てみます．

n=5以降で比較的変化が緩やかになっていることがわかります．ただ，n=5以降の値の変化に注目すると面白い現象が見えました．

グラフから分かるように，nの増加に伴って誤差が単調に減少するわけではなく，増減が見られます．また，nが奇数の時に比較的誤差の小さい結果が得られていることもわかります．これは，ガウス点を求めるsubroutineの仕様である，nが奇数のときには強制的に0をガウス点の座標として追加するという効果が表れているものだと考えられます．