プログラミングと数値計算

ガウス=クロンロッド求積法

コメントする

ガウス=クロンロッド求積法と呼ばれる数値積分法です。

  1. ガウス=クロンロッド求積法の説明
  2. fortran90によるプログラム
    1. 15次ガウス=クロンロッド求積(\([-1\sim 1]\))
    2. 15次ガウス=クロンロッド求積(\([a\sim b]\))
    3. 61次ガウス=クロンロッド求積(\([a\sim b]\))
  3. 参考文献

説明

ガウス=ルジャンドル求積法の補助的な存在、というとらえ方でいいと思います。
通常はガウス=ルジャンドル求積法の誤差推定のために使われるらしいです。

積分
\(
\displaystyle \int_{-1}^{1} f(x) dx
\)
を、
\(
\displaystyle \int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{2N+1} \omega_i f(x_i)+\mathcal{O}(x^{3N+1})
\)
として近似します。ここで、\(\omega_i\)は点\(x_i\)での重みを意味します。
\(\omega_i, x_i\)は任意に決められず、\(N\)をいくつにするかによって強制的に決定されます。

ガウス=クロンロッド求積法は次数(\(2N+1\))を用いたとき、被積分関数が
\(x^{3N+1}\)以下の多項式であれば厳密な積分値を返しますCHAPTER 5 Numerical Quadrature。
Nは次数。分点数と一致。

重要な事ですが、被積分関数が多項式ではない場合、あまりよくない積分結果になってしまいます。
例えば、積分
\(
\displaystyle \int_0^b \sqrt{x} dx
\)
を考えると、被積分関数が多項式で書けないので、精度はガクッと落ちます。
この理由は、被積分関数の微分が\(x=0\)で発散するためであり、高次の積分法になるほど顕著に表れます。

また、ルンゲ現象に代表される、高次で発生する様々な問題を排除するため、通常は積分区間を複数に分けて、低次の積分法を組み合わせて計算を行います。

[adsense1]

fortran90によるプログラム

以下のプログラムは15次ガウス=クロンロッド求積法のプログラムです。
これは、積分
\(
\displaystyle \int_{-1}^{1}e^{-x} dx = 2\sinh(1) \approx 2.350402387287602913764764
\)
を計算するプログラムとなっています。

program main
  implicit none
  integer::i,N
  double precision,allocatable::x(:),w(:)
  double precision::f,s
  external::f
 
  N=15
  allocate(x(1:N),w(1:N)); x=0d0; w=0d0
  call GaussKronrod15(x,w)
  s=0d0
  do i=1,N
     s=s+w(i)*f(x(i))
  enddo  
  write(6,*)s

  deallocate(x,w)  
 
  stop
end program main

function f(x)
  implicit none
  double precision,intent(in)::x
  double precision::f
 
  f=exp(-x)

  return
end function f

subroutine GaussKronrod15(x,w)
  !Gauss-Kronrod Quadrature Nodes and Weights
  !http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
  implicit none
  double precision,intent(out)::x(1:15),w(1:15)
 
  integer::i
  integer,parameter::N=15
 
  x=0d0; w=0d0
 
  x( 8) = 0d0
  x( 9) = 2.077849550078984676006894037732449d-1
  x(10) = 4.058451513773971669066064120769615d-1
  x(11) = 5.860872354676911302941448382587296d-1
  x(12) = 7.415311855993944398638647732807884d-1
  x(13) = 8.648644233597690727897127886409262d-1
  x(14) = 9.491079123427585245261896840478513d-1
  x(15) = 9.914553711208126392068546975263285d-1  
  do i=1,7
     x(i)=-x(N-i+1)
  enddo
 
  w( 8) = 2.094821410847278280129991748917143d-1
  w( 9) = 2.044329400752988924141619992346491d-1
  w(10) = 1.903505780647854099132564024210137d-1
  w(11) = 1.690047266392679028265834265985503d-1
  w(12) = 1.406532597155259187451895905102379d-1
  w(13) = 1.047900103222501838398763225415180d-1
  w(14) = 6.309209262997855329070066318920429d-2
  w(15) = 2.293532201052922496373200805896959d-2
  do i=1,7
     w(i)=w(N-i+1)
  enddo

  return
end subroutine GaussKronrod15

実行結果は
2.35040238728760
となり、厳密値
2.350402387287602913
と全桁一致します。

区間変更

定義そのままでは\([-1\sim 1]\)なので、区間を変更しましょう。
区間変更は\(y=\frac{b-a}{2}x+\frac{b+a}{2}\)と変数変換をすれば、

\(
\begin{align}
\int_a^b f(x) dx &= \int_{-1}^{1} f(y) \frac{b-a}{2}dy \\
&\approx \sum_{i=1}^N \left(\frac{b-a}{2}\omega_i\right) f(\frac{b-a}{2}x+\frac{b+a}{2}) \\
&= \sum_{i=1}^N \omega_i’ f(x_i’)
\end{align}
\)
と掛けます。これをプログラムします。
ここでは、積分
\(
\displaystyle \int_{1}^{2.4}e^{-x} dx = e^{-1}-e^{-2.4} \approx 0.27716148788202981
\)
を計算します。

program main
  implicit none
  integer::i,N
  double precision,allocatable::x(:),w(:)
  double precision::f,s
  external::f
 
  N=15
  allocate(x(1:N),w(1:N)); x=0d0; w=0d0
  call GaussKronrod15ab(1d0,2.4d0,x,w)
  s=0d0
  do i=1,N
     s=s+w(i)*f(x(i))
  enddo  
  write(6,*)s

  deallocate(x,w)  
 
  stop
end program main

function f(x)
  implicit none
  double precision,intent(in)::x
  double precision::f
 
  f=exp(-x)

  return
end function f

subroutine GaussKronrod15ab(a,b,x,w)
  !Gauss-Kronrod Quadrature Nodes and Weights
  !http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
  implicit none
  double precision,intent(in)::a,b
  double precision,intent(out)::x(1:15),w(1:15)
 
  integer::i
  integer,parameter::N=15
 
  x=0d0; w=0d0
 
  x( 8) = 0d0
  x( 9) = 2.077849550078984676006894037732449d-1
  x(10) = 4.058451513773971669066064120769615d-1
  x(11) = 5.860872354676911302941448382587296d-1
  x(12) = 7.415311855993944398638647732807884d-1
  x(13) = 8.648644233597690727897127886409262d-1
  x(14) = 9.491079123427585245261896840478513d-1
  x(15) = 9.914553711208126392068546975263285d-1  
  do i=1,7
     x(i)=-x(N-i+1)
  enddo
 
  w( 8) = 2.094821410847278280129991748917143d-1
  w( 9) = 2.044329400752988924141619992346491d-1
  w(10) = 1.903505780647854099132564024210137d-1
  w(11) = 1.690047266392679028265834265985503d-1
  w(12) = 1.406532597155259187451895905102379d-1
  w(13) = 1.047900103222501838398763225415180d-1
  w(14) = 6.309209262997855329070066318920429d-2
  w(15) = 2.293532201052922496373200805896959d-2
  do i=1,7
     w(i)=w(N-i+1)
  enddo
 
  x=0.5d0*((b-a)*x+(a+b))
  w=0.5d0*(b-a)*w
 
  return
end subroutine GaussKronrod15ab

となります。

実行結果は
0.277161487882030
となり、厳密値
0.27716148788202981
と全桁一致します。数値計算なので、最後のところは四捨五入されています。

61次ガウス=クロンロッド求積法

61次ガウス=クロンロッド求積法のコードも置いておきます。
これほどの高次は余り使わず、低次のものを組み合わせたほうが良いことを再び書いておきます。

積分
\(
\displaystyle \int_{-3}^{20}e^{-x} dx = e^{3}-e^{-20} \approx 20.0855369211265141
\)
を計算します。

program main
  implicit none
  integer::i,N
  double precision,allocatable::x(:),w(:)
  double precision::f,s
  external::f
 
  N=61
  allocate(x(1:N),w(1:N)); x=0d0; w=0d0
  call GaussKronrod61ab(-3d0,20d0,x,w)
  s=0d0
  do i=1,N
     s=s+w(i)*f(x(i))
  enddo  
  write(6,*)s

  deallocate(x,w)  
 
  stop
end program main

function f(x)
  implicit none
  double precision,intent(in)::x
  double precision::f
 
  f=exp(-x)

  return
end function f

subroutine GaussKronrod61ab(a,b,x,w)
  !Gauss-Kronrod Quadrature Nodes and Weights
  !http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
  implicit none
  double precision,intent(in)::a,b
  double precision,intent(out)::x(1:61),w(1:61)
 
  integer::i
  integer,parameter::N=61
 
  x=0d0; w=0d0
 
  x(31:61) = &
       (/0.000000000000000000000000000000000d0  &
       , 5.147184255531769583302521316672257d-2 &
       , 1.028069379667370301470967513180006d-1 &
       , 1.538699136085835469637946727432559d-1 &
       , 2.045251166823098914389576710020247d-1 &
       , 2.546369261678898464398051298178051d-1 &
       , 3.040732022736250773726771071992566d-1 &
       , 3.527047255308781134710372070893739d-1 &
       , 4.004012548303943925354762115426606d-1 &
       , 4.470337695380891767806099003228540d-1 &
       , 4.924804678617785749936930612077088d-1 &
       , 5.366241481420198992641697933110728d-1 &
       , 5.793452358263616917560249321725405d-1 &
       , 6.205261829892428611404775564311893d-1 &
       , 6.600610641266269613700536681492708d-1 &
       , 6.978504947933157969322923880266401d-1 &
       , 7.337900624532268047261711313695276d-1 &
       , 7.677774321048261949179773409745031d-1 &
       , 7.997278358218390830136689423226832d-1 &
       , 8.295657623827683974428981197325019d-1 &
       , 8.572052335460610989586585106589439d-1 &
       , 8.825605357920526815431164625302256d-1 &
       , 9.055733076999077985465225589259583d-1 &
       , 9.262000474292743258793242770804740d-1 &
       , 9.443744447485599794158313240374391d-1 &
       , 9.600218649683075122168710255817977d-1 &
       , 9.731163225011262683746938684237069d-1 &
       , 9.836681232797472099700325816056628d-1 &
       , 9.916309968704045948586283661094857d-1 &
       , 9.968934840746495402716300509186953d-1 &
       , 9.994844100504906375713258957058108d-1 /)

  do i=1,30
     x(i)=-x(N-i+1)
  enddo
  w(31:61)= &
       (/5.149472942945156755834043364709931d-2 &
       , 5.142612853745902593386287921578126d-2 &
       , 5.122154784925877217065628260494421d-2 &
       , 5.088179589874960649229747304980469d-2 &
       , 5.040592140278234684089308565358503d-2 &
       , 4.979568342707420635781156937994233d-2 &
       , 4.905543455502977888752816536723817d-2 &
       , 4.818586175708712914077949229830459d-2 &
       , 4.718554656929915394526147818109949d-2 &
       , 4.605923827100698811627173555937358d-2 &
       , 4.481480013316266319235555161672324d-2 &
       , 4.345253970135606931683172811707326d-2 &
       , 4.196981021516424614714754128596976d-2 &
       , 4.037453895153595911199527975246811d-2 &
       , 3.867894562472759295034865153228105d-2 &
       , 3.688236465182122922391106561713597d-2 &
       , 3.497933802806002413749967073146788d-2 &
       , 3.298144705748372603181419101685393d-2 &
       , 3.090725756238776247288425294309227d-2 &
       , 2.875404876504129284397878535433421d-2 &
       , 2.650995488233310161060170933507541d-2 &
       , 2.419116207808060136568637072523203d-2 &
       , 2.182803582160919229716748573833899d-2 &
       , 1.941414119394238117340895105012846d-2 &
       , 1.692088918905327262757228942032209d-2 &
       , 1.436972950704580481245143244358001d-2 &
       , 1.182301525349634174223289885325059d-2 &
       , 9.273279659517763428441146892024360d-3 &
       , 6.630703915931292173319826369750168d-3 &
       , 3.890461127099884051267201844515503d-3 &
       , 1.389013698677007624551591226759700d-3/)
  do i=1,30
     w(i)=w(N-i+1)
  enddo
 
  x=0.5d0*((b-a)*x+(a+b))
  w=0.5d0*(b-a)*w
 
  return
end subroutine GaussKronrod61ab

実行結果
20.0855369211265
厳密値
20.0855369211265141

[adsense2]

fortran90で、任意の次数の求積点を知りたければ、f90のプログラムが
KRONROD Gauss-Kronrod Quadrature Rules -Source Codes in
Fortran90
にて公開されています。

参考文献

Gauss-Kronrod Quadrature Nodes and Weights

関連記事

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です