Clenshaw–Curtis求積法

クレンショウ-カーティス(Clenshaw–Curtis)型の数値積分公式です。

[2]の論文”IS GAUSS QUADRATURE BETTER THAN CLENSHAW–CURTIS?”の結論だけかいつまみますと、

・クレンショウ-カーティスの求積法は、基本的にはガウス-ルジャンドル求積法と同じ性能
・クレンショウ-カーティスの求積法は実装が簡単

ということです。
本気で計算精度を求めようとするならばガウス-ルジャンドル求積法、簡単に、でもかなりの精度で求めたければクレンショウ-カーティスの求積法なのだと思います。
ガウス-ルジャンドル求積法は分点の位置と重みの計算が簡単ではないのです。

積分
\(
\displaystyle \int_{-1}^1 f(x) dx \sim \sum_{k=0}^{n} \omega_k f(x_k)
\)

Nが偶数の時、分点\(x_k\)と重み\(\omega_k\)は
\(
\displaystyle x_k=\cos\left(\frac{k\pi}{N}\right),\;\;\;k=0,1,\cdots,N
\)

\(
\begin{align}
\omega_0 &=\omega_N =\frac{1}{N^2 -1} \\
\omega_s &=\omega_{N-s} = \frac{4}{n^2}\sum^{n/2}_{j=0}{\prime\prime }\frac{1}{1-4j^2}\cos\left(\frac{2j\pi s}{n}\right),\;\;\;s=0,1,\cdots,\frac{N}{2}
\end{align}
\)

ここで、\(\sum{\prime\prime }\)は和の最初と最後の項に\(\frac{1}{2}\)を乗じることを意味します。

N点の分点に対し、誤差のオーダーは\(O(N+1)\)です。
\(x_k\)はチェビシェフ多項式\(T_n(x)\)の極値点になっています。

利点


・重みの計算に掛かる時間が\(O(N log N)\).
・2N点での積分値を求める際にN点で計算した評価点を用いることが出来る。

ガウス求積法と比べて


・ガウス型の重みの計算に掛かる時間は\(O(N^2)\).ただし、分点と重みは一度求めてしまえばずっと使いまわせるので、さほど重要な問題ではありません。
・ガウス-ルジャンドル求積法の分点を\(N\)点利用し、それに\(N\)点追加することでガウス-クロンロッド求積法が使えます。ガウス-ルジャンドルはN点の求積で\(O(x^{2N})\)、ガウス-クロンロッドは2N点の求積で\(O(x^{3N})\)の多項式まで正確に表すことが出来ます。

結論

本気で計算精度が欲しいならばガウス-ルジャンドル求積法。
簡単に、でもかなりの精度で求めたければクレンショウ-カーティスの求積法なのだと思います。

プログラム


最適化は特にしていませんが、fortran90ではこんな感じになります。

program main
  implicit none
  integer::i,j,N
  double precision::a,b,s
  double precision,allocatable::x(:),w(:)
  double precision,parameter::pi=acos(-1d0)
  double precision::f
  external::f
 
  a=0d0
  b=pi
   
  N=11
  allocate(x(1:N),w(1:N))
  call clenshaw_curtis_xw(N,x,w)
  w=w*(b-a)*0.5d0
  x=(x*(b-a)+(b+a))*0.5d0

  s=0d0
  do i=1,N
     s=s+f(x(i))*w(i)
  enddo
  write(6,*)s

  stop
end program main
   
function f(x)
  implicit none
  double precision::f,x
  double precision::pi=dacos(-1d0)
 
  f=sin(x)*x
 
  return
end function

subroutine clenshaw_curtis_xw(N,x,w)
  implicit none
  integer,intent(in)::N
  double precision,intent(out)::x(1:*),w(1:*)
  double precision,parameter::pi=3.14159265358979323846264338327950288d0
  integer::i,j,M,k
  double precision::c
 
  if(mod(N,2).eq.0)then
     write(6,*)"N must be odd"
     stop
  endif
 
  M=N-1
 
  do i=0,M
     x(i+1)=cos(i*pi/dble(M))
  enddo
  w(1)=1d0/((M-1d0)*(M+1d0))
  w(M+1)=w(1)

  do i=1,M/2
     k=i+1
     w(k)=0.5d0*(1d0+cos(pi*i)/((1d0-M)*(1d0+M)))
     c=2d0*pi*i/M
     do j=1,M/2-1
        w(k)=w(k)+cos(c*j)/((1d0-2d0*j)*(1d0+2d0*j))
     enddo
     w(k)=w(k)*4d0/M
     
     w(M-i+1)=w(k)
  enddo

  return
end subroutine clenshaw_curtis_xw

参考文献

[1]P. J. Davis, P. Rabinowitz著、森正武訳 『計算機による数値積分法』p. 131-132,初版 (1980)
[2]LLOYD N. TREFETHEN, IS GAUSS QUADRATURE BETTER THAN CLENSHAW–CURTIS?


「Clenshaw–Curtis求積法」への5件のフィードバック

  1. コメント失礼します.

    分点XkはNが偶数のときしか計算できないのですか?

    1. 少なくとも本稿で紹介しているやり方では出来ない、そういう意味です。
      奇数の時も出来るかもしれませんが、私は調べてはいません。

  2. 「分点と重みは一度求めてしまえばずっと使いまわせるので、さほど重要な問題ではありません」とあるのですが,つまり,まず低い精度の数値積分で用いた分点の値を,さらに高い精度の数値積分に用いることができるということですか?

  3. 例えば,
    N=4のとき
    分点は,0, 0.707, 6.125E-17, -0.707, -1
    N=6のとき
    分点は,0, 0.866, 0.5, 6.125E-17, 0-0.5, -0.866, -1
    となりますが

    つまり,0,-1,6.125E-17の値が使いまわせるとう意味でしょうか?
    低次のときの全点が使い回せるわけではないですよね?

    1. ゆうさん

      コメントに気付くのが遅くなりすみません。

      >「分点と重みは一度求めてしまえばずっと使いまわせるので、さほど重要な問題ではありません」とあるのですが,つまり,まず低い精度の数値積分で用いた分点の値を,さらに高い精度の数値積分に用いることができるということですか?

      上記の「」は、本文中ではガウス求積法に関する話です。Clenshaw–Curtis求積法の話ではありません。
      Clenshaw–Curtis求積法は高い精度に使いまわすことは出来ませんが,
      ガウス-ルジャンドル求積法は分点や重みを追加するだけで, ガウス-クロンロッド求積法という別の次数の計算結果が得られます。

      同じ分点の数で複数の積分を行いたい場合を想定した説明であり、1個の積分をするために毎回分点と重みを計算する必要はない、だから「さほど重要な問題では」ないということです。

ゆう へ返信する コメントをキャンセル

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