最速・高精度の数値積分

無償で公開されている1次元数値積分コードで優秀なルーチンは何か?を見つけます。

何よりも早さが足りない人へ。
言語はFortranです。

まとめ


QUADPACK(ルーチン名dqag, key=2,6)
http://www.netlib.org/quadpack/
が広い場合で使えて優秀です。しかもパブリックドメイン!
直ぐ使えるように私が変更を加えたQUADPACKプログラムが下の方にあります。

もう少し問題に対して最適化をしたければ以下のチャートに従ってください。

もくじ

  1. 比較対象について(QUADPACK, NUMPAC, Ooura’s MSP)
  2. ライセンスについて
  3. 比較方法
  4. 結果
  5. QUADPACKプログラム
  6. ベンチマーク時のプログラム
  7. 参考文献

比較対象(QUADPACK, NUMPAC, Ooura’s MSP)


どんな関数が入ってきても計算時間と精度を両立させる万能な数値積分手法は存在しません。一長一短があります。

この一長一短は非常に極端です。
被積分関数に特定の性質があるのであればその方法を使うべきです。
特定の積分範囲、被積分関数が(重み関数)×(多項式)で表現されるのであれば、ガウス求積法が最速で高精度です。
被積分関数が端点で特異性(例えば端点で\(x^{-1/2}\))を持つのであれば二重指数関数型数値積分が最速で高精度です。

ここでは、無償公開されている自動積分コードのみを比較します。
使っている手法が実行したい積分とマッチすれば自動積分は最善の方法です。

比較するコードは以下のものです。

比較コード
パッケージ
or 作成者(ルーチン名)
積分戦略 積分の手法
NUMPAC(aqnn9d)
http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/
適応型 ニュートン・コーツ9点則
NUMPAC(rombgd)
http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/
非適応型 ロンバーグ積分
Ooura(intcc)
http://www.kurims.kyoto-u.ac.jp/~ooura/intcc-j.html
非適応型 クレンショウ-カーティス積分
Ooura(intde)
http://www.kurims.kyoto-u.ac.jp/~ooura/index-j.html
非適応型 二重指数関数型数値積分
QUADPACK(dqag,key=2)
http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html
適応型 7-15点ガウス-ルジャンドル、
ガウス-クロンロッド求積法
QUADPACK(dqag,key=6)
http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html
適応型 30-61点ガウス-ルジャンドル、
ガウス-クロンロッド求積法
Sikino(romberg)
http://slpr.sakura.ne.jp/qp/nc-integral/
非適応型 ロンバーグ積分
Sikino(de)
http://slpr.sakura.ne.jp/qp/de-fourmula/
非適応型 二重指数関数型数値積分
Sikino(rkf451)
http://slpr.sakura.ne.jp/qp/runge-kutta-ex/
適応型 ルンゲ-クッタ-フェールベルグ法による積分

 適応型というのは、被積分関数の形に依って分点を増やす減らすを決める自動積分法です。
非適応型というのは、被積分関数の形に依らず分点を増やす減らすを決める自動積分法です。


スポンサーリンク


ライセンスについて


各プログラムの利用条件に関して注記しておきます。
私の理解であるので、解釈に間違いがあるかもしれません。最終的な判断は個人の判断にお任せします

NetlibにあるQUADPACK
http://www.netlib.org/quadpack/
パブリックドメイン。Fortran数値積分用のパッケージ。
著作権表示無しに、商用非商用問わず複製、利用、改変、配布して良いもの。
ただし、所有権や人格権などを侵してはならない。

私がNetlibにあるQUADPACKをfortran90用に改変したものを置いておきます。下のQUADPACKプログラムへ。

people.sc.fsu.eduにあるQUADPACK
GNU LGPLライセンス。Fortran数値積分用のパッケージ。

! Licensing:
!
! This code is distributed under the GNU LGPL license.

とあり、http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html
にコードがあります。こちらの方が使うまでが簡単ですが、ライセンスがあることが面倒。

NUMPAC

NUMPACを利用した論文には,NUMPAC利用の記述を記載してください。
また,著者から修正プログラムの提供があったときには置き換えてください。
NUMPAC邦文マニュアル – XML

とあります。著作権表示や複製、利用、改変、配布に関する記述は無いので、これらをするべきではありません。

Ooura’s Mathematical Software Packages
プログラムのreadmeの中に

copyright
Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
You may use, copy, modify this code for any purpose and
without fee. You may distribute this ORIGINAL package.

とあります。翻訳すれば、
”著作権はTakuya OOURAにあります。どんな目的に対しても無償でコードを使用、複製、改変して良いです。そしてオリジナルパッケージを配布して良いです。”
とあります。
例えば上の説明にあてはまらないため、許諾なしに改変したものを配布してはいけません

本サイト(シキノート)のプログラム
私のは最終的な出版物や何かに、私の著作権表示、利用した旨を記述してもらえれば商用非商用であっても複製、利用、改変、配布して良いです。
論文に使用した場合、表示は要りません。
改変したものを配布したい場合、著作権は私にある、しかし責任は改変を行った者にあることを注記してください。改変したものについてもそれによって生じた責任は負いません。

すなわち本当に自由に使えるのはNetlibで公開されているQUADPACKだけです。

比較方法


Kahaner’s 21 test problemsより、問題を抜粋して比較します。
Kahaner’s 21 test problemsとは21個のベンチマーク用の問題で、全て解析解が存在するものです。
これが解けたら大した数値積分プログラムだ!ということです。
被積分関数に特異性や発散などを含んでいる問題です。

より、以下の10問を選んで比較します。問題の抜粋は私が独断で選びました。

結果


計算結果を並べます。比較は、

  • 積分プログラムに積分区間と精度を要求し、その精度に達したか?
  • 要求精度を実現するまでに被積分関数が評価された回数

を比較します。また、適応型の積分手法に対して、実際にどの点が評価されているかを知りたかったため、
実際に評価された点達に対し矩形波で畳み込みを行い、最大値を1にするように規格化し比較しました。
これが左から1枚目の画像です。
2枚目は、真値との相対誤差を要求精度\(\varepsilon\)の関数として、
3枚目は、要求精度が達成された場合の被積分関数の評価回数です。
3枚目が一番重要な情報です。

この問題(2),(3)は被積分関数に不連続性、特異性がある場合の問です。
(2)の有効な手法は適応型の数値積分法であることが分かります。シンプルな問題ですが、大域的に精度を決める自動積分では3-4桁程度にとどまり、全く計算できていません。
(3)は左の端点で一回微分が発散する関数です。最も有効な手法は変数変換型の数値積分、二重指数関数型数値積分法で、次いで適応型の積分法です。一回微分が発散が意味するのは、テーラー展開による誤差評価で、誤差項が発散することを意味しています。


(3)はラグランジュ補間が難しい関数です。ラグランジュ補間が難しい場合、チェビシェフ多項式による補間が有効です。
(9)は振動型の関数です。sin(x)のテーラー展開が難しいように、これもチェビシェフ多項式による補間が有効です。


(13)は早く振動する関数です。関数に節の数が多いため次数の高い手法が必要、もしくは分割するしかありません。
(14)は振動型の関数です。sin(x)のテーラー展開が難しいように、これもチェビシェフ多項式による補間が有効です。


(17)はx=0回りだけに値が存在します。適応型でなければ、局在部分を評価するために要らない部分を増やさなければならない無駄が生じてしまいます。
(18)は程々の振動型の関数です。チェビシェフ多項式による補間が有効です。


(20)はローレンツ型の関数です。ローレンツ型はラグランジュ補間ではうまく補間されません。
(21)は突如値が大きくなる地点が存在する関数です。うまく変化する点を見つけないと的外れな計算結果になります。

ちなみに


二重指数関数型数値積分が有効な例として
\(
\displaystyle \int_0^{1}\sin(1/\sqrt{x})/\sqrt{x}
\)

という例題が挙げられます。難しい事には変わりませんが二重指数関数型数値積分だけが出来るのではないことをここに注記しておきます。

二重指数関数型数値積分の凄いところは、適当に作った私のプログラムでさえ有名な積分プログラムと張り合えるという点でしょう。

非適応型の数値積分法の重みは被積分関数の形に依存しないで決まることを確認しておきましょう。いつでもこんな感じの重みになります。


まとめです。

問題番号 被積分関数の特徴 最も有効な積分法 次いで有効な積分法
(2) 不連続性 NUMPAC(aqnn9d) QUADPACK(dqag, key=2)
(3) 端点特異性 Ooura(intde) NUMPAC(aqnn9d)
(5) ローレンツ Ooura(intcc) QUADPACK(dqag, key=6)
(9) 振動関数 Ooura(intcc) QUADPACK(dqag, key=6)
(13) 高振動 Ooura(intcc) QUADPACK(dqag, key=6)
(14) 局在 NUMPAC(aqnn9d) Ooura(intcc)
(17) 局在、振動 Ooura(intcc) QUADPACK(dqag, key=6)
(18) 振動 Ooura(intcc) QUADPACK(dqag, key=6)
(20) ローレンツ Ooura(intcc) QUADPACK(dqag, key=6)
(21) 突如変化する関数 NUMPAC(aqnn9d) NUMPAC(rombgd)

ある精度で積分の問題を解きたい時に、最高のパフォーマンスを出したければ以下のチャートに従うと良いと思います。ライセンスの問題を含めて、QUADPACKが優秀でしょう。

QUADPACKのkey=2がよいか、key=6がよいかは良く分かりません。上の表にはkey=6が優秀に見えますが、\(\sqrt{x}\)の積分などでは圧倒的にkey=2が優秀だからです。

QUADPACKプログラム


私が手を加えているので、オリジナルと比較し相違点が生じるかもしれません。
また、使用するにあたっては、私の作成したものと同様、元にしたものが私の物である、と表示さえしていただければ、商用非商用問わず複製、利用、改変、配布してくださって結構です。ただし、責任は負いません。
問題ないと思いますが、もし公開しているプログラムに著作権的にダメな点があればお知らせください。判明次第、ページを非公開にし、早急に対応いたします。

Quadpackのコード
http://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag.f90

プログラムのサンプルコード

program main
  implicit none
  integer::ier
  double precision::a,b,s,eps,g
  external::g
 
  a=0d0
  b=1d0
  eps=1d-12

  call dqag_k2(g,a,b,eps,s,ier)
  write(6,*)s,ier
  ! ier = 0 : success
  ! ier > 0 : fail

  call dqag_k6(g,a,b,eps,s,ier)
  write(6,*)s,ier
 
  stop
end program main

function g(x)
  implicit none
  double precision::g,x
 
  g=sqrt(x)

  return
end function g

g  : 被積分関数\(g(x)\)
a,b : 積分区間\([a,b]\)
eps : 要求する相対誤差
s  : \(g(x)\)の数値積分値
ier : ier=0 成功、収束した。
     ier>0 失敗、収束しなかった。

本来のquadpackには色々な引数がありますが、それを省略するためにdqag_k2、dqag_k6というサブルーチンを挟んでいます。これはquadpack_dqag.f90の中に入っているので参照してください。
dqag_k2は7-15点ガウス-ルジャンドル、ガウス-クロンロッド求積法を用いた適応型の数値積分法です。
dqag_k6は30-61点ガウス-ルジャンドル、ガウス-クロンロッド求積法を用いた適応型の数値積分法です。

コンパイルと実行は

$ gfortran quadpack_dqag.f90 main.f90
$ ./a.out
  0.66666666666666718                0
  0.66666666666666796                0
$

でokです。積分値が0の時、相対誤差の条件は満たされないためier=2が返ってきます。
問題に応じて無視するかどうか決めてください。

もしも値が存在するはずなのにier=1が帰ってくるようであれば、サブルーチンdqag_k2、dqag_k6の中のlimitの値を増やしてください。こんなふうに。

  ier=1
  epsabs=-1d0  
  !limit=200
  limit=1000
  lenw=limit*4

2017/11/02 追記)
複素平面上の点\(a=(a_x,a_y)\)から、\(b=(b_x,b_y)\)に至る経路を直線で結ぶ、複素関数\(g(z)\)の線積分は以下のコードです。

ただし、多価関数は不可能です。偏角の範囲はfortran90では\([-\pi\lt \theta \le \pi]\)出なければなりません。

Quadpackのコード
http://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90

プログラムのサンプルコード

program main
  implicit none
  integer::ier
  double precision::eps
  complex(kind(0d0))::a,b,s
  complex(kind(0d0)),external::g
 
  a=dcmplx(0d0,3d0)
  b=dcmplx(1d0,2d0)
  eps=1d-12

  call cqag_k2(g,a,b,eps,s,ier)
  write(6,*)s,ier
 
  stop
end program main

function g(z)
  implicit none
  complex(kind(0d0))::g
  complex(kind(0d0)),intent(in)::z

  g=sqrt(z)
 
  return
end function g
スポンサーリンク

ベンチマーク用プログラム


自分の持っているプログラムのベンチマークを行いたい場合のプログラムを公開します。

参考文献


QUADPACK Netlibhttp://www.netlib.org/quadpack/
QUADPACK_DOUBLE Numerical Integration http://people.sc.fsu.edu/~jburkardt/f_src/quadpack_double/quadpack_double.html
NUMPAC http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/
Ooura’s Mathematical Software Packages http://www.kurims.kyoto-u.ac.jp/~ooura/i
シキノート
http://slpr.sakura.ne.jp/qp/nc-integral/

kahaner’s 21 test problemに関して
Takemitsu Hasegawa, Susumu Hibino, Yohsuke Hosoda and Ichizo Ninomiya,
“An extended doubly-adaptive quadrature method based on the combination of the Ninomiya and the FLR schemes”

幸谷智紀、鈴木千里 “補外法に基づく数値積分法の実装と性能評価”


コメントを残す

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