適応型シンプソン法

適応シンプソン法は適応シンプソン則とも呼ばれ、 1962 年に GF Kuncir が提案した数値積分の手法です。^[1]これはおそらく印刷物に登場した最初の再帰的適応型数値積分アルゴリズムですが、^{[2]現在では}、ガウス–クロンロッド積分法やクレンショウ–カーティス積分法に基づくより現代的な適応型手法が一般的に好まれています。適応シンプソン法は、シンプソン則を使用して定積分を計算することで得られる誤差の推定値を使用します。誤差がユーザー指定の許容値を超える場合、このアルゴリズムでは積分区間を 2 つに分割し、各区間に適応シンプソン法を再帰的に適用します。この手法は、関数が3 次関数で十分に近似される場所で関数評価の回数が少なくなるため、通常、合成シンプソン則よりもはるかに効率的です。

シンプソン則は、積分関数が3次以下の多項式である場合に正確な補間求積法です。リチャードソン外挿法を用いると、6つの関数値に対するより正確なシンプソン推定値と、3つの関数値に対するより精度の低い推定値に、補正を適用することで組み合わせることができます。したがって、得られた推定値は5次以下の多項式に対して正確です。 $S(a,m)+S(m,b)$ $S(a,b)$ $[S(a,m)+S(m,b)-S(a,b)]/15$

数学的手順

用語の定義

^{JN Lyness [3]}によって提案された、区間の細分化をいつ停止するかを決定する基準は、 $|S(a,m)+S(m,b)-S(a,b)|<15\varepsilon$

ここで、は中点を持つ区間であり、シンプソンの定理によって与えられる、、はそれぞれ、、の推定値であり、は区間の望ましい最大許容誤差です。 $[a,b]$ $m={\frac {a+b}{2}}$ $S(a,b)\,\!$ $S(a,m)\,\!$ $S(m,b)$ ${\textstyle \int _{a}^{b}f(x)\,dx}$ ${\textstyle \int _{a}^{m}f(x)\,dx}$ ${\textstyle \int _{m}^{b}f(x)\,dx}$ $\varepsilon$

注記、。 $\varepsilon _{i+1}=\varepsilon _{i}/2$

手順

適応型シンプソン法を実行するには、次の操作を行います。の場合、積分を近似するために使用されるシンプソン則の和にとを加算し、それ以外の場合はの代わりにとを使用して同じ操作を実行します $|S(a,m)+S(m,b)-S(a,b)|<15\varepsilon _{i}$ $S(a,m)$ $S(m,b)$ ${\textstyle \left|S{\left(a,{\frac {m-a}{2}}\right)}+S{\left({\frac {m-a}{2}},m\right)}-S(a,m)\right|<15\varepsilon _{i+1}}$ ${\textstyle \left|S{\left(m,{\frac {b-m}{2}}\right)}+S{\left({\frac {b-m}{2}},b\right)}-S(m,b)\right|<15\varepsilon _{i+1}}$ $|S(a,m)+S(m,b)-S(a,b)|<15\varepsilon _{i}$

数値的考察

適応型シンプソン法では、入力によっては収束が遅く、許容値がアンダーフローして無限ループが発生します。この問題を防ぐ簡単な方法としては、深さ制限を追加する（CサンプルやMcKeemanのように）、浮動小数点演算で $ε /2\neq ε$ を検証する、またはその両方（Kuncirのように）などがあります。区間サイズはローカルマシンのイプシロンに近づき、 $a = bとなる場合もあります$

ライネスの1969年の論文には、この問題をより具体的に扱う「修正4」が含まれている: ^[3]^{: 490–2}

初期区間を $[A, B]$ とし、元の許容値を $ε 0$ とする。
$各部分区間[a, b]$ について、誤差推定値 $D (a, b)をと定義する$ 。E $= 180$ $ε$ $0$ $/ ($ $B$ $-$ $A$ $)$ $と$ 定義する。この場合、元の終了基準は $D$ $\leq$ $E$ となる。 ${\textstyle {\frac {12}{a-b}}[S^{(2)}(a,b)-S(a,b)]={\frac {1}{4}}[f(a)-4f(lm)+6f(m)-4f(rm)+f(b)]}$
D ( a , m ) ≥ D ( a , b )の場合、四捨五入レベルに達したか、f ⁽⁴⁾が区間内にゼロ点を持つかのいずれかである。許容誤差ε _0をε′ ₀に変更する必要がある。
- 再帰ルーチンは、現在の区間のDレベルを返す必要があります。ルーチン静的変数 $E' = 180 ε' 0 / (B - A)が定義され、$ Eに初期化されます。
- （修正4 i、ii）区間でさらに再帰が使用される場合：
  1. 四捨五入されたと思われる場合は、E'を $D (a, m)$ に変更します。^[a]
  2. それ以外の場合は、E'を $max(E 、 D (a 、 m))$ に調整します。
- 調整にはある程度のコントロールが必要です。許容範囲の大幅な増加やわずかな減少は抑制する必要があります。
区間全体にわたって有効ε′ _{0を計算するには：}
- E 'が変化する各 $x i を$ $、$ $($ $x$ $i$ $,ε$ $i$ $'$ $)$ のペアの配列に対数で記録する。最初の要素は $($ $a$ $,$ $ε'0$ $)$ となる。
- 実際のε _effは、区間の幅に応じて重み付けされたすべてのε′ ₀の算術平均です。
ある区間の電流E'がEよりも高い場合、5次の加速/補正は適用されない：^[b]
- 終了基準の「15」要素は無効になっています。
- 訂正用語は使用しないでください。

イプシロン増加法は、ルーチンを「ベストエフォート」モードで使用できるようにします。初期許容値がゼロの場合、ルーチンは最も正確な答えを取得し、実際のエラーレベルを返します。^[3]^{: 492}

サンプルコードの実装

以下に示す一般的な実装手法は、区間 $[$ $a$ $,$ $b$ $]とともに$ $f(a), f(b), f(m)を下位に渡すことです。親レベルで$ $S$ $($ $a$ $,$ $b$ $)$ を評価するために使用されるこれらの値は、部分区間にも再び使用されます。これにより、各再帰呼び出しのコストが入力関数の評価回数6回から2回に削減されます。使用されるスタック領域のサイズは、再帰の層に比例します

Python

以下はPythonでの適応型シンプソン法の実装です

from __future__ import division # python 2 互換# Racket から翻訳された「構造化された」適応型バージョンdef _quad_simpsons_mem ( f , a , fa , b , fb ): """シンプソンの定理を評価し、再利用するために m と f(m) を返します""" m = ( a + b ) / 2 fm = f ( m ) return ( m , fm , abs ( b - a ) / 6 * ( fa + 4 * fm + fb ))    

     
    
          
      
                   

def _quad_asr ( f , a , fa , b , fb , eps , whole , m , fm ): """    適応型シンプソン則の効率的な再帰実装。    区間の開始、中間、終了時の関数値が保持されます。    """ lm , flm , left = _quad_simpsons_mem ( f , a , fa , m , fm ) rm , frm , right = _quad_simpsons_mem ( f , m , fm , b , fb ) delta = left + right - whole if abs ( delta ) <= 15 * eps : return left + right + delta / 15 return _quad_asr ( f , a , fa , m , fm , eps / 2 , left , lm 、flm ) + \
            _quad_asr ( f 、m 、fm 、b 、fb 、eps / 2 、右、rm 、frm )         
    



             
            
          
         
               
                       

def quad_asr ( f , a , b , eps ): """適応型シンプソン則を使用して、最大誤差 eps で a から b まで f を積分します。""" fa , fb = f ( a ), f ( b ) m , fm , whole = _quad_simpsons_mem ( f , a , fa , b , fb ) return _quad_asr ( f , a , fa , b , fb , eps , whole , m , fm )    
    
        
            
             

math import sin print ( quad_asr ( sin , 0 , 1 , 1e-09 ))から

C

これは、fと求積法の計算の冗長な評価を回避する、C99での適応型シンプソン法の実装です。数値問題に対する3つの「単純な」防御策すべてが含まれています

#include <math.h>   // fabs と sin のインクルードファイル#include <stdio.h> // printf と perror のインクルードファイル#include <errno.h> 
 
 

/** 適応型シンプソン則、再帰コア */ 
float adaptiveSimpsonsAux ( float ( * f )( float ), float a , float b , float eps , float whole , float fa , float fb , float fm , int rec ) { float m = ( a + b ) / 2 , h = ( b - a ) / 2 ; float lm = ( a + m ) / 2 , rm = ( m + b ) / 2 ; // 深刻な数値的問題: 収束しませんif (( eps / 2 == eps ) || ( a == lm )) { errno = EDOM ; return whole ; } float flm = f ( lm ), frm = f ( rm ); float left = ( h / 6 ) * ( fa + 4 * flm + fm );フロート右= ( h / 6 ) * ( fm + 4 * frm + fb );フロートデルタ=左+右-全体;        
                                    
                   
                 
    
                  
               
              
             
           

    if ( rec <= 0 && errno != EDOM ) errno = ERANGE ; // 深度制限が浅すぎます// Lyness 1969 + Richardson 外挿。記事を参照してくださいif ( rec <= 0 || fabs ( delta ) <= 15 * eps ) return left + right + ( delta ) / 15 ; return adaptiveSimpsonsAux ( f , a , m , eps / 2 , left , fa , fm , flm , rec -1 ) + adaptiveSimpsonsAux ( f , m , b , eps / 2 , right , fm , fb , frm , rec -1 ); }            
    
           
             
               
                   


/** 適応型シンプソン則ラッパー
* (キャッシュされた関数評価を埋めます) */ 
float adaptiveSimpsons ( float ( * f )( float ), // 積分する関数ポインタfloat a , float b , // 区間 [a,b] float epsilon , // 誤差許容度int maxRecDepth ) { // 再帰キャップerrno = 0 ; float h = b - a ; if ( h == 0 ) return 0 ; float fa = f ( a ), fb = f ( b ), fm = f (( a + b ) / 2 ); float S = ( h / 6 ) * ( fa + 4 * fm + fb ); return adaptiveSimpsonsAux ( f , a , b , epsilon , S , fa , fb , fm , maxRecDepth ); }       
                                
                                 
                              
      
         
         
               
           
             


/** 使用例 */ 
#include <stdlib.h> // 敵対的な例 (rand 関数) の場合static int callcnt = 0 ; static float sinfc ( float x ) { callcnt ++ ; return sinf ( x ); } static float frand48c ( float x ) { callcnt ++ ; return drand48 (); } int main () { // I を 0 から 2 までの sin(x) の積分としますfloat I = adaptiveSimpsons ( sinfc , 0 , 2 , 1e-5 , 3 ); printf ( "integrate(sinf, 0, 2) = %lf \n " , I ); // 結果を出力しますperror ( "adaptiveSimpsons" ); // 成功しましたか? (depth=1 は浅すぎます) printf ( "(%d 評価) \n " , callcnt ); 
    
         
         
  
    
           
        
                       
     

    callcnt = 0 ; srand48 ( 0 ); I = adaptiveSimpsons ( frand48c , 0 , 0.25 , 1e-5 , 25 ); // ランダム関数printf ( "integrate(frand48, 0, 0.25) = %lf \n " , I ); perror ( "adaptiveSimpsons" ); // 収束しませんprintf ( "(%d 評価) \n " , callcnt ); return 0 ; }

この実装は、オークリッジ国立研究所^[4]をはじめとするプロジェクトにおけるX線レーザーシミュレーション用のC++レイトレーサーに組み込まれています。ORNL版は、呼び出しカウンタ、様々なデータ型のテンプレート、多次元積分用のラッパーなどが追加され、拡張されています。^[4]

ラケット

これは、 Racketにおける適応型シンプソン法の実装であり、動作ソフトウェアコントラクトを用いて実装されている。エクスポートされた関数は、与えられた関数fの不定積分を計算する。^[5]

;; ----------------------------------------------------------------------------- 
;; インターフェース、コントラクト付き
( provide/contract 
[adaptive-simpson ( ->i (( f ( -> real? real? )) ( L real? ) ( R ( L ) ( and/c real? ( >/c L )))) ( #:epsilon ( ε real? )) ( r real? )) ] )  
               
         
        

;; ----------------------------------------------------------------------------- 
;; 実装
( define ( adaptive-simpson f L R #:epsilon [ε . 000000001] ) ( define f@L ( f L )) ( define f@R ( f R )) ( define-values ( M f@M whole ) ( simpson-1call-to-f f L f@L R f@R )) ( asr f L f@L R f@R ε whole M f@M ))       
     
     
           
           

;; 「効率的な」実装
( define ( asr f L f@L R f@R ε whole M f@M ) ( define-values ( leftM f@leftM left* ) ( simpson-1call-to-f f L f@L M f@M )) ( define-values ( rightM f@rightM right* ) ( simpson-1call-to-f f M f@M R f@R )) ( define delta* ( - ( + left* right* ) whole )) ( cond [ ( <= ( abs delta* ) ( * 15 ε )) ( + left* right* ( / delta* 15 )) ] [else ( define epsilon1 ( / ε 2 )) ( + ( asr f L f@L M f@M epsilon1 left* leftM f@leftM ) ( asr f M f@M R f@R epsilon1 right* rightM f@rightM )) ] ))          
              
           
        
  
               
         
                       
                      

;; 区間の半分を評価する (1 func eval) 
( define ( simpson-1call-to-f f L f@L R f@R ) ( define M ( mid L R )) ( define f@M ( f M )) ( values M f@M ( * ( / ( abs ( - R L )) 6 ) ( + f@L ( * 4 f@M ) f@R ))))      
      
     
                 

( ( mid L R ) ( / ( + L R ) 2. )を定義します。)

注釈

^ オリジナルの4iではE'を上げることしか言及されていません。しかし、後のテキストでは、E'を大きく下げることができると述べられています
^ これは、単純なケースでは、許容値/間隔幅のアンダーフローにも当てはまる可能性があります。

参考文献

^ ab GF Kuncir (1962)、「アルゴリズム103：シンプソンのルール積分器」、Communications of the ACM、5（6）：347、doi：10.1145 / 367766.368179
^ abより ODEソルバーを彷彿とさせる、初期の非再帰適応積分器については、S. Henriksson (1961)「Contribution no. 2: Simpson numeric integration with variable length of step」、BIT Numerical Mathematics、1 :290を参照。
^ abcdef JN Lyness (1969)、「適応型シンプソン求積ルーチンに関する注記」、Journal of the ACM、16 (3): 483– 495、doi : 10.1145/321526.321537
^ ab ベリル、マーク A (2016 年 9 月 27 日)。「RayTrace ミニアプリ: src/AtomicModel/interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9」。ORNL GitLab。
^ Felleisen, Matthias. 「[racket] adaptive simpson integration」. Racketメーリングリスト「ユーザー」 . 2018年9月26日閲覧。
^ マッキーマン、ウィリアム・マーシャル（1962年12月1日）「アルゴリズム145：シンプソン則による適応型数値積分」Communications of the ACM . 5 (12): 604. doi : 10.1145/355580.369102 .

外部リンク

適応型シンプソン則のモジュール

[4] オリジナルの4iではE'を上げることしか言及されていません。しかし、後のテキストでは、E'を大きく下げることができると述べられています

[5] これは、単純なケースでは、許容値/間隔幅のアンダーフローにも当てはまる可能性があります。

[kuncir-1] GF Kuncir (1962)、「アルゴリズム103：シンプソンのルール積分器」、Communications of the ACM、5（6）：347、doi：10.1145 / 367766.368179

[henriksson-2] より ODEソルバーを彷彿とさせる、初期の非再帰適応積分器については、S. Henriksson (1961)「Contribution no. 2: Simpson numeric integration with variable length of step」、BIT Numerical Mathematics、1 :290を参照。

[lyness-3] JN Lyness (1969)、「適応型シンプソン求積ルーチンに関する注記」、Journal of the ACM、16 (3): 483– 495、doi : 10.1145/321526.321537

[ornl-6] ベリル、マーク A (2016 年 9 月 27 日)。「RayTrace ミニアプリ: src/AtomicModel/interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9」。ORNL GitLab。

[7] Felleisen, Matthias. 「[racket] adaptive simpson integration」. Racketメーリングリスト「ユーザー」 . 2018年9月26日閲覧。

[8] マッキーマン、ウィリアム・マーシャル（1962年12月1日）「アルゴリズム145：シンプソン則による適応型数値積分」Communications of the ACM . 5 (12): 604. doi : 10.1145/355580.369102 .