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?

奥多摩探索(工場と廃墟への行き方)

奥多摩は東京にありながら自然豊かです。

ここでは、奥多摩にある工場(奥多摩工業㈱)と廃墟(奥多摩湖ロープウェイ、みとうさんぐち駅)への道順を記述していきます。

奥多摩工業㈱ 氷川工場

概要


奥多摩工業㈱ 氷川工場は石灰の採掘、販売を行う稼働中の工場です。
稼働中なのですが、一部壁が剥がれ落ちていたり、増築に増築を重ねたようなSFを漂わす感じが非常に魅力です。
(奥多摩工業 -wikipedia)

場所



綺麗に見える場所は日原川から望む展望でしょう。
上の地図上の、緑色の場所から望めます。
こんな感じ

工場脇を流れる日原川は澄み切っています。

また、工場の内部に立ち入ることはできませんが、工場を横切るように一般道が存在しています。
この一般道は地図上の赤いポイントから赤いポイントです。

上流から望む氷川工場の写真がwikipediaにありました。この場所は工場のすぐ上にある釣り堀から撮ったと思われます。

untitled

[adsense1]

奥多摩湖ロープウェイ、みとうさんぐち駅

概要


奥多摩湖ロープウェイは小河内観光開発株式会社が運行を行っていましたが、1975年に正式に運行休止され、現在では事実上の廃墟になっています。
ロープウェイ自体は1962年-1966年まで操業されていました。
現在、地元自治体の申し出によって立ち入りは禁止されているようです(wikipedia情報 小河内観光開発)。

50年近く経った現在でもコンクリート作りの建物であるため、しっかりと残っています。

場所


場所は奥多摩周遊道路、川野駐車場のすぐ横の階段を昇った先にあります。

(2017年5月確認)日祝の場合
・バス停「奥多摩駅」発 → バス停「陣屋」(じんや)降車(乗車時間30分)
発車時刻 2番バスターミナル、奥12 (7:25、10:05、13:35、16:35)
  コ : 留浦経由小菅の湯行 奥12
  小 : 大菩薩峠東口経由小菅の湯行 奥12
  奥多摩駅、バス時刻表 http://transfer.navitime.biz/bus-navi/pc/diagram/BusDiagram?orvCode=00042103&course=0000431622&stopNo=1

・バス停「陣屋」発 → 奥多摩駅行き(乗車時間30分)
発車時刻 奥12(8:56、11:46、15:07、18:07)

陣屋、バス時刻表http://bus.ekitan.com/rosen/Rp700?t=0&b=310615&f=0&b2=310558&com=31&cn=%90%BC%93%8C%8B%9E%83o%83X

※後から気が付きましたが、バスは1区間遠い『深山橋』で降りる/乗ると奥12だけでなく、奥09,奥10,奥11でも良いようです。この場合、一時間に一本バスが出ているようです。

奥多摩駅時刻表
http://transfer.navitime.biz/bus-navi/pc/diagram/BusDiagram?orvCode=00042103&course=0000431622&stopNo=1

深山橋→奥多摩駅行きの時刻表
http://transfer.navitime.biz/bus-navi/pc/diagram/BusDiagram?orvCode=00042635&course=0000430401&stopNo=19&orvName=%E6%B7%B1%E5%B1%B1%E6%A9%8B&poleId=

[adsense2]