最速・高精度の数値積分

無償で公開されている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”

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

フリーソフトの紹介

全てwindows10上で動きます。
Microsoft Powerpoint, Word以外は無料です。

  1. キャプチャ -Setuna2
  2. ランチャー -RocketDuck
  3. メール管理 -Thuderbird
  4. テキストエディタ -Sublime Text
  5. 文書作成 -Microsoft Word, -LaTeX
  6. グラフ作成 -gnuplot
  7. PDFビューワ- SumatraPDF
  8. 画像圧縮 -Caesium portable
  9. 画像編集 -Microsoft Powerpoint, -GIMP
  10. GIF画像作成、編集 -Giam
  11. 画像ビューワ -Leeyes
  12. 音楽再生 -VLCメディアプレイヤー
  13. 音楽ダウンロード -Download helper(firefoxアドオン)
  14. 動画再生 -GOMplayer
  15. 動画変換、音楽の抽出 -Craving Explorer
  16. mp3の音量調節 -MP3Gain
  17. 動画キャプチャ -oCam
  18. 動画編集 -Aviutl
  19. PCの情報取得 -Speccy
  20. PC間のデータ移動 -IP Messenger
  21. PCのデータ削除 -Bleachbit portable
  22. PCのウイルス除去 -Norton Power Eraser
  23. アプリケーション情報 -Process Explorer
  24. 広告除去 -Adblock Plus
  25. 仮想環境 -Vmware Player
  26. オフラインゲーム -Track Mania Forever

1. キャプチャ


Setuna2
http://www.clearunit.com/clearup/setuna2/

説明

画面上の任意の範囲を画像として保存するソフトウェアです。
画像保存だけでなく、キャプチャした領域は半透明で常に画面最前列に現れます。
そのため、アプリケーション間を行ったり来たりすることなく欲しい情報を表示させておけます。
現在サイトが閉鎖されている…ようです?


2. ランチャー


RocketDuck
https://rocketdock.com/

説明

ランチャーです。Macにある感じのやつです。画像の上に出てる枠を追加します。
良く使う実行ファイルやフォルダへのショートカットをデスクトップ上以外に表示させることが出来ます。


3. メール管理


Thuderbird
https://www.mozilla.org/ja/thunderbird/

説明

複数のメールアドレス(例えばyahoo,gmail,hotmailなど他社サービスであっても)をまとめて管理できます。
インターネットブラウザを立ち上げないでメールを見られます。


4. テキストエディタ


Sublime Text
https://www.sublimetext.com/3

説明

プログラムを書くためのエディタです。


5. 文書作成


Microsoft Word, LaTeX

説明

文章を書くためのソフトウェアです。
数式や記号のようなものを表現したければLaTeXが良いです。
LaTeXは導入までが難しいです。導入はEasyTeXが比較的簡単でしょう。
EasyTeX


6. グラフ作成


gnuplot
http://www.gnuplot.info/

説明

グラフを描く為のソフトウェアです。
windows用に作られたwxtターミナルはアンチエイリアスが聞いているので、滑らかなグラフが表示されます。
eps,gif画像も作ることが出来ます。


7. PDFビューワ


SumatraPDF
https://www.sumatrapdfreader.org/free-pdf-reader.html

説明

PDFファイルを開くためのソフトウェアです。
PDFファイルを見るだけしか出来ませんが、起動が高速で、動作中に固まったことはありません。
PDFビューワとして有名なソフトウェアAdobe Readerがありますが、起動までに数秒~数十秒かかり動作も重いです。


スポンサーリンク



8. 画像圧縮


Caesium portable
https://saerasoft.com/caesium/

説明

jpg,png等の画像を圧縮し,容量を減少させます。
インストールなしに使うことができ、手軽に使えます。


9. 画像編集


Microsoft Powerpoint
GIMP
https://www.gimp.org/

説明

・Powerpoint
本来プレゼンテーションを作るためのソフトです。
画像や写真に文字を加えたい時などに、ppt上で編集→図として保存を行うと比較的簡単に編集できます。

・GIMP
画像編集用ソフトです。
フォトショップ並みの編集が可能ですが、慣れ、知識、技術が必要ですので、初心者向けではありません。


10. GIF画像作成、編集


Giam
http://furumizo.net/tsu/giamd.htm

説明

GIF画像の編集ソフトです。
GIF画像を作ったり、編集を加えたり、コマ送りとして見ることが出来るソフトです。
ある特定の場面を画像として保存することも出来ます。


11. 画像ビューワ


Leeyes
http://www3.tokai.or.jp/boxes/leeyes/

説明

画像、または圧縮ファイルにおさめられた画像を解凍することなく見ることが出来る画像ビューワです。
プラグインを入れることによりrarなどの様々な画像形式に対応することが出来ます。


12. 音楽再生


VLCメディアプレイヤー
https://www.videolan.org/vlc/index.ja.html

説明

音楽、動画再生ソフトです。


13. 音楽,動画ダウンロード


Download helper(firefoxアドオン)
https://addons.mozilla.org/ja/firefox/addon/video-downloadhelper/

説明

firefoxのアドオンで、あらゆるサイトの動画のダウンロードが可能です。
3つの風船がトレードマークで、動画配信サイトに行きこの風船をクリックし、動画ダウンロードを行います。
動画のダウンロードは、著作権法に抵触しないようにしましょう。
※ブラウザに付随するコンソールでも、メディアファイルのダウンロードが可能です。


14. 動画再生


GOMplayer
http://www.gomplayer.jp/

説明

動画再生ソフトです。
コマ送り、スロー再生、逆再生、左右反転、上下反転が出来ます。
広告が邪魔ですが、これは環境設定→一般→スキンからクラシック[フォルダ]を選ぶと広告は消えます。


15. 動画変換、音楽の抽出


Craving Explorer
https://www.crav-ing.com/

説明

本来はインターネットブラウザです。
Craving explorerに付属する、オフライン状況下でも使える動画変換機能がシンプルで優秀です。
以前は動画、音楽ダウンロードを売りにしていたブラウザでしたが、著作権の問題からダウンロード機能が除去されていたようです。
2017年現在、ホームページを訪れると再び動画をダウンロードできる仕様になったようです。


16. mp3の音量調節


MP3Gain
http://mp3gain.sourceforge.net/

説明

mp3ファイルの音量を変えるソフトです。
複数ファイルの音量の均一化をすることも出来ます。
エンコードするわけではないので、操作は一瞬で完了します。


17. 動画キャプチャ


oCam
http://ohsoft.net/eng/

説明

画面を動画としてキャプチャします。
録画後に広告が出ますが、録画した動画に広告は出てきません。


18. 動画編集


Aviutl
http://spring-fragrance.mints.ne.jp/aviutl/

説明

avi動画の編集を行えます。
全く別の2つの動画から比較動画を作ったり、文字を入力したりできます。
初心者向けではありません。編集には慣れが必要です。
プラグインを用いれば、mp4を読み込ませ、編集し、mp4で出力させることが出来ます。
詳しい説明は
VIPで初心者がゲーム実況するには https://www18.atwiki.jp/live2ch/pages/165.html

ニコニコ動画まとめwiki http://nicowiki.com/AviUtl.html
をご覧ください。

スポンサーリンク


19. PCの情報取得


Speccy
https://www.piriform.com/speccy

説明

PCの構成情報を綺麗なウインドウで教えてくれるだけのソフトウェアです。
例えば、OS, CPU, メモリ(DDR3など)、モニターの型番、マウスの情報、IPアドレス等々です。
スマートです。


20. PC間のデータ移動


IP Messenger
https://ipmsg.org/

説明

ルーターを介してLANでつながったPC間のデータ移動を行います。
インターネット環境は必要ありません。ネット接続もしません。
高速、セキュリティともに良いと思います。開発から20年経過した現在(2017/09/08)でも更新がなされています。


21. PCのデータ削除


Bleachbit portable
https://www.bleachbit.org/

説明

PCに存在するゴミファイルを消去して空き容量を確保します。
ブラウザのキャッシュ、一時ファイル、データを揃えて空き容量を確保したりできます。
アプリケーションごとにどういったファイルを消すか指定することが出来ます。


22. PCのウイルス除去


Norton Power Eraser
https://security.symantec.com/nbrt/npe.aspx?lcid=1041

説明

感染したPCからウイルスを除去する強力なツールです。
徹底的にデータを消すため、無害なソフトを巻き込むこともあるようです。
システムの復元でも消えない、hao123の除去も出来ます。


23. アプリケーション情報


Process Explorer
https://docs.microsoft.com/en-us/sysinternals/downloads/process-explorer

説明

不審な動きをするアプリケーションを特定するために使えます。


24. 広告除去


Adblock Plus
https://addons.mozilla.org/ja/firefox/addon/adblock-plus/

Element Hiding helper for Adblock Plus
https://addons.mozilla.org/ja/firefox/addon/elemhidehelper/

説明

インターネットのサイトに表示させる広告を除去することが出来ます。
スマホ(Android)用のAdblockブラウザもありますが、動作が遅いです。
firefox、Chromeにはプラグインがあります。

Adblock Plusは画像広告に強く、
Element Hiding helper for Adblock Plusはテキスト広告に強いです。


25. 仮想環境


Vmware Player
https://www.vmware.com/products/player/playerpro-evaluation.html

説明

windows上にlinux環境を構築できます。
あくまでwindows上の仮想環境なので、linuxが壊れてしまった場合やlinuxを消したい場合も単にwindows上からフォルダを消すだけokです。


26. オフラインゲーム


Track Mania Forever
http://trackmaniaforever.com/
無料版ダウンロード http://www.4gamer.net/games/048/G004822/20080423014/

説明

オフラインで遊べるレーシングゲームです。
あるコースを時間内に走り抜けることが目的のゲームです。
コースエディタがあるので、自分で作ることが可能です。