Adaptacyjna metoda Simpsona

Adaptacyjna metoda Simpsona , zwana także adaptacyjną regułą Simpsona , to metoda całkowania numerycznego zaproponowana przez GF Kuncira w 1962 r. Jest to prawdopodobnie pierwszy rekurencyjny algorytm adaptacyjny całkowania numerycznego, który pojawił się w druku, chociaż bardziej nowoczesne metody adaptacyjne oparte na kwadraturze Gaussa – Kronroda i kwadratura Clenshawa – Curtisa są obecnie ogólnie preferowane. Metoda Adaptive Simpsona wykorzystuje oszacowanie błędu, jaki otrzymujemy przy obliczaniu całki oznaczonej za pomocą reguły Simpsona . Jeśli błąd przekracza tolerancję określoną przez użytkownika, algorytm wzywa do podzielenia przedziału całkowania na dwie części i zastosowania adaptacyjnej metody Simpsona do każdego podprzedziału w sposób rekurencyjny. Technika ta jest zwykle znacznie bardziej wydajna niż złożona reguła Simpsona, ponieważ wykorzystuje mniej ocen funkcji w miejscach, w których funkcja jest dobrze aproksymowana przez funkcję sześcienną .

Reguła Simpsona jest interpolacyjną regułą kwadraturową, która jest dokładna, gdy całka jest wielomianem stopnia trzeciego lub niższego. Korzystając z ekstrapolacji Richardsona $) + S (m, b)}$ dokładniejsze oszacowanie Simpsona dla sześciu wartości funkcji jest połączone z mniej dokładnym oszacowaniem ${\ Displaystyle S (a, b)}$ dla trzech wartości funkcji przez zastosowanie poprawki ${\ styl wyświetlania [S(a,m)+S(m,b)-S(a,b)]/15}$ . Uzyskane oszacowanie jest więc dokładne dla wielomianów stopnia piątego lub mniejszego.

Procedura matematyczna

Definiowanie warunków

Kryterium określania, kiedy przestać dzielić przedział, sugerowane przez JN Lynessa, jest

$${\ Displaystyle | S (a, m) + S (m, b) -S (a, b) | <15 \ varepsilon}$$

gdzie ${\ Displaystyle [a, b]}$ jest przedziałem z punktem środkowym ${\ Displaystyle m = {\ Frac {ba} {2}}}$ , podczas gdy ${\ Displaystyle S (a, b) \, \!}$ , ${\ Displaystyle S (a, m) \, \!}$ i ${\ Displaystyle S (m, b)}$ podane przez regułę Simpsona to oszacowania ${\ textstyle \ int _ {a} ^ {b} f (x) \, dx$ , ${\ textstyle \ int _ {a} ^ {m} f (x) \, dx} i$ ∫ $\ textstyle \ int _ {m} ^ {b} f (x) \, dx}$ odpowiednio, a $jest$ maksymalną tolerancją błędu dla przedziału.

Uwaga, ${\ Displaystyle \ varepsilon _ {i + 1} = \ varepsilon _ {i}/2}$ .

Kroki proceduralne

Aby wykonać adaptacyjną metodę Simpsona, wykonaj następujące czynności: if ${\ Displaystyle | S (a, m) + S (m, b) -S (a, b) | <15 \ varepsilon _ {i}} dodaj$ S $\ Displaystyle S (a, m)}$$,$ reguł które są używane do przybliżenia całki, w przeciwnym razie wykonaj tę samą ${\ textstyle \ lewo | S {\ lewo (a, {\ Frac {ma} {2}} \ prawo)} + S {\ lewo ({\ Frac {ma} {2}}, m\right)}-S(a,m)\right|<15\varepsilon _{i+1}}$ i ${\textstyle \ lewo|S {\ lewo (m, {\ Frac {bm} {2}} \ prawo)} + S {\ lewo ({\ Frac {bm} {2}}, b\right)}-S(m,b)\right|<15\varepsilon _{i+1}}$ zamiast ${\ Displaystyle | S (a, m) + S (m, b) -S (a, b) | <15 \ varepsilon _ {i}}$ .

Rozważanie liczbowe

Niektóre dane wejściowe nie zostaną szybko zbieżne w adaptacyjnej metodzie Simpsona, co spowoduje niedopełnienie tolerancji i utworzenie nieskończonej pętli. Proste metody ochrony przed tym problemem obejmują dodanie ograniczenia głębokości (jak w próbce C i McKeeman), sprawdzenie, czy ε /2 ≠ ε w arytmetyce zmiennoprzecinkowej lub jedno i drugie (jak Kuncir). Rozmiar interwału może również zbliżyć się do epsilon maszyny lokalnej , dając a = b .

Artykuł Lychee z 1969 roku zawiera „Modyfikację 4”, która odnosi się do tego problemu w bardziej konkretny sposób:

• Niech początkowy przedział będzie [ A , B ] . Niech pierwotną tolerancją będzie ε ₀ .
• Dla każdego podprzedziału [ a , b ] zdefiniuj D ( a , b ) , oszacowanie błędu jako ${\ textstyle {\ Frac {12} {ab}} [S ^ {(2)} ( a,b)-S(a,b)]={\frac {1}{4}}[f(a)-4f(lm)+6f(m)-4f(rm)+f(b)]}$ . Zdefiniujmy ₀ E = 180 ε / ( B - A ) . Pierwotne kryteria zakończenia miałyby wówczas postać D ≤ E .
• ₀ Jeśli D ( a , m ) ≥ D ( a , b ) , albo osiągnięto poziom zaokrąglenia, albo w przedziale znaleziono zero dla f ^{(4) .} Konieczna jest zmiana tolerancji ε na ε′ . ₀
  • Podprogramy rekurencyjne muszą teraz zwracać poziom D dla bieżącego interwału. Zmienna rutynowo-statyczna ₀ E' = 180 ε' / ( B - A ) jest zdefiniowana i zainicjowana jako E .
  • (Modyfikacja 4 i, ii) Jeśli dalsza rekurencja jest używana w przedziale:
    1. Jeśli wydaje się, że zaokrąglenie zostało osiągnięte, zmień E' na D ( a , m ) .
    2. W przeciwnym razie dostosuj E' do max( E , D ( a , m )) .
  • Konieczna jest pewna kontrola regulacji. Należy powstrzymać znaczne wzrosty i niewielkie spadki tolerancji.
• ₀ Aby obliczyć efektywną ε′ w całym przedziale:
  • Zarejestruj każde x _i , przy którym E' zmienia się w tablicę par ( x _i , ε _i ' ) . Pierwszym wpisem powinno być ₀ ( a , ε′ ) .
  • ₀ Rzeczywisty ε _eff jest średnią arytmetyczną wszystkich ε′ ważoną szerokością przedziałów.
• Jeżeli bieżące E' dla przedziału jest większe niż E , to przyspieszenie/poprawka piątego rzędu nie miałaby zastosowania:
  • Współczynnik „15” w kryteriach zakończenia jest wyłączony.
  • Nie należy stosować terminu korygującego.

Manewr podnoszenia epsilon pozwala na użycie procedury w trybie „najlepszego wysiłku”: przy zerowej początkowej tolerancji, procedura spróbuje uzyskać najbardziej precyzyjną odpowiedź i zwróci rzeczywisty poziom błędu.

Przykładowe implementacje kodu

Powszechną techniką implementacji pokazaną poniżej jest przekazywanie w dół f( a ), f( b ), f( m ) wraz z przedziałem [ a , b ] . Te wartości, użyte do oceny S ( a , b ) na poziomie nadrzędnym, zostaną ponownie użyte dla podprzedziałów. W ten sposób zmniejsza się koszt każdego wywołania rekurencyjnego z 6 do 2 ocen funkcji wejściowej. Rozmiar używanej przestrzeni stosu pozostaje liniowy w stosunku do warstwy rekurencji.

Pyton

Oto implementacja adaptacyjnej metody Simpsona w Pythonie .

    

     
    
          
      
                   

         
    



             
            
          
         
               
                       

    
    
        
            
             

   
 0   from  __future__  import  Division  # python 2 compat  # "structured" adaptacyjna wersja, przetłumaczona z Racket  def  _quad_simpsons_mem  (  f  ,  a  ,  fa  ,  b  ,  fb  ):  """Ocenia regułę Simpsona, zwracając również m i f(m) do ponownego użycia """  m  =  (  za  +  b  )  /  2  fm  =  fa  (  m  )  powrót  (  m  ,  fm  ,  abs  (  b  -  za  )  /  6  *  (  fa  +  4  *  fm  +  fb  ))  def  _quad_asr  (  fa  ,  za  ,  fa  ,  b  ,  fb  ,  eps  ,  whole  ,  m  ,  fm  ):  """  Efektywna rekurencyjna implementacja adaptacyjnej reguły Simpsona.  Zachowane są wartości funkcji na początku, w środku, na końcu przedziału.  """  lm  ,  flm  ,  left  =  _quad_simpsons_mem  (  f  ,  a  ,  fa  ,  m  ,  fm  )  rm  ,  frm  ,  right  =  _quad_simpsons_mem  (  f  ,  m  ,  fm  ,  b  ,  fb  )  delta  =  lewa  +  prawa  -  całość  jeśli  abs  (  delta  )  <=  15  *  eps  :  powrót  w lewo  +  prawo  +  delta  /  15  powrót  _quad_asr  (  f  ,  a  ,  fa  ,  m  ,  fm  ,  eps  /  2  ,  lewo  ,  lm  ,  flm  )  +  \  _quad_asr  (  f  ,  m  ,  fm  ,  b  ,  fb  ,  eps  /  2  ,  prawo  ,  rm  ,  frm  )  def  quad_asr  (  f  ,  a  ,  b  ,  eps  ):  """Całkuj f od a do b za pomocą adaptacyjnej reguły Simpsona z maksymalnym błędem eps."""  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  ,  całość  ,  m  ,  fm  )  z  matematyki  import  sin  print  (  quad_asr  (  sin  ,  ,  1  ,  1e-09  )) 

C

Oto implementacja adaptacyjnej metody Simpsona w C99, która pozwala uniknąć zbędnych ocen f i obliczeń kwadraturowych. Obejmuje wszystkie trzy „proste” sposoby obrony przed problemami numerycznymi.

 
 
 


        
                                    
                   
                 
    
                  
               
              
             
           

       0         
    
       0    
             
               
                   




       
                                
                                 
                              
      0
         
       0  0
               
           
             



 
    0
         
         
  
    
        0   
        
                       
     

      0 0
       0    
     
                            
     
     0
 #include  <math.h>  // plik dołączania dla fabs i sin  #include  <stdio.h>  // plik dołączany dla printf i perror  #include  <errno.h>  /** Adaptacyjna reguła Simpsona, rdzeń rekurencyjny */  float  adaptiveSimpsonsAux  (  float  (  *  f  )(  float  ),  float  a  ,  float  b  ,  float  eps  ,  float  całość  ,  float  fa  ,  float  fb  ,  float  fm  ,  int  rec  )  {  float  m  =  (  a  +  b  )  /  2  ,  h  =  (  b  -  a  )  /  2  ;  float  lm  =  (  za  +  m  )  /  2  ,  rm  =  (  m  +  b  )  /  2  ;  // poważny problem liczbowy: nie będzie zbieżny  if  ((  eps  /  2  ==  eps  )  ||  (  a  ==  lm  ))  {  errno  =  EDOM  ;  zwrócić  w całości  ;  }  float  flm  =  fa  (  lm  ),  frm  =  fa  (  rm  );  pływak  w lewo  =  (  h  /  6  )  *  (  fa  +  4  *  flm  +  fm  );  pływak  w prawo  =  (  h  /  6  )  *  (  fm  +  4  *  frm  +  fb  );  float  delta  =  lewa  +  prawa  -  całość  ;  if  (  rec  <=  &&  errno  !=  EDOM  )  errno  =  ERANGE  ;  // granica głębokości zbyt płytka  // Lyness 1969 + ekstrapolacja Richardsona; zobacz artykuł   if  (  rec  <=  ||  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  ) ;  }  /** Adaptacyjne opakowanie reguły Simpsona  * (wypełnia oceny funkcji w pamięci podręcznej) */  float  adaptiveSimpsons  (  float  (  *  f  )(  float  ),  // funkcja ptr do całkowania  float  a  ,  float  b  ,  // interwał [a,b]  float  epsilon  ,  // tolerancja błędu  int  maxRecDepth  )  {  // limit rekurencji  errno  =  ;  liczba zmiennoprzecinkowa  h  =  b  -  a  ;  jeśli  (  h  ==  )  powrót  ;  float  fa  =  fa  (  za  ),  fb  =  fa  (  b  ),  fm  =  fa  ((  a  +  b  )  /  2  );  float  S  =  (  godz  /  6  )  *  (  fa  +  4  *  fm  +  fb  );  return  adaptiveSimpsonsAux  (  f  ,  a  ,  b  ,  epsilon  ,  S  ,  fa  ,  fb  ,  fm  ,  maxRecDepth  );  }  /** przykład użycia */  #include  <stdlib.h>  // dla wrogiego przykładu (funkcja Rand)  static  int  callcnt  =  ;  static  float  sinfc  (  float  x  )  {  callcnt  ++  ;  zwróć  sinf  (  x  );  }  static  float  frand48c  (  float  x  )  {  callcnt  ++  ;  zwróć  drand48  ();  }  int  main  ()  {  // Niech I będzie całką sin(x) od 0 do 2  float  I  =  adaptiveSimpsons  (  sinfc  ,  ,  2  ,  1e-5  ,  3  );  printf  (  "integrate(sinf, 0, 2) = %lf  \n  "  ,  I  );  // wypisz wynik  perror  (  "adaptiveSimpsons"  );  // Czy to się udało? (głębokość = 1 jest zbyt płytka)   printf  (  "(%d ocen)  \n  "  ,  callcnt  );  wywołanie  =  ;  srand48  (  );  I  =  adaptacyjne Simpsony  (  frand48c  ,  ,  0.25  ,  1e-5  ,  25  );  // losowa funkcja  printf  (  "integrate(frand48, 0, 0.25) = %lf  \n  "  ,  I  );  perror  (  "adaptacyjne Simpsony"  );  // nie będzie zbieżny  printf  (  "(%d ocen)  \n  "  ,  callcnt  );  powrót  ;  } 

Ta implementacja została włączona między innymi do śledzenia promieni C++ przeznaczonego do symulacji lasera rentgenowskiego w Oak Ridge National Laboratory . Wersja ORNL została wzbogacona o licznik połączeń, szablony dla różnych typów danych i opakowania do integracji w wielu wymiarach.

Rakieta

Oto implementacja adaptacyjnej metody Simpsona w Racket z behawioralnym kontraktem na oprogramowanie. Wyeksportowana funkcja oblicza całkę nieoznaczoną dla danej funkcji f .

  
               
         
        



       
     
     
           
           


          
              
           
        
  
               
         
                       
                      


      
      
     
                 

         ;; -------------------------------------------------- ---------------------------   ;; interfejs, z kontraktem   (  dostarcz/umowa  [adaptive-simpson  (  ->i  ((  f  (  ->  real?  real?  ))  (  L  real?  )  (  R  (  L  )  (  and/c  real?  (  >/c  L  ) )))  (  #:epsilon  (  ε  real?  ))  (  r  real?  ))  ]  )  ;; -------------------------------------------------- ---------------------------   ;; implementacja   (  zdefiniuj  (  adaptive-simpson  f  L  R  #:epsilon  [ε  .  000000001]  )  (  zdefiniuj  f@L  (  f  L  ))  (  zdefiniuj  f@R  (  f  R  ))  (  zdefiniuj-wartości  (  M  f@M  całość  )  (  simpson-1call-to-f  f  L  f@L  R  f@R  ))  (  asr  f  L  f@L  R  f@R  ε  całe  M  f@M  ))  ;; „wydajna” implementacja   (  zdefiniuj  (  asr  f  L  f@L  R  f@R  ε  całe  M  f@M  )  (  zdefiniuj wartości  (  leftM  f@leftM  left*  )  (  simpson-1call-to-f  f  L  f@L  M  f@M  ))  (  zdefiniuj wartości  (  prawyM  f@prawyM  prawy*  )  (  simpson-1call-to-f  f  M  f@M  R  f@R  ))  (  zdefiniuj  delta*  (  -  (  +  lewy*  prawy*  )  całość  ))  (  warunek  [  (  <=  (  abs  delta*  )  (  *  15  ε  ))  (  +  lewy*  prawy*  (  /  delta*  15  ))  ]  [else  (  zdefiniuj  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  ))  ]  ))  ;; oceń połowę przedziału (1 func eval)   (  zdefiniuj  (  simpson-1call-to-f  f  L  f@L  R  f@R  )  (  zdefiniuj  M  (  mid  L  R  ))  (  zdefiniuj  f@M  (  f  M  ))  (  wartości  M  f@M  (  *  (  /  (  abs  (  -  R  L  ))  6  )  (  +  f@ L  (  *  4  f@M  )  f@R  ))))  (  zdefiniuj  (  mid  L  R  )  (  /  (  +  L  R  )  2.  )) 

Powiązane algorytmy

• Henriksson (1961) to nierekurencyjny wariant reguły Simpsona. „Dostosowuje się”, całkując od lewej do prawej i dostosowując szerokość interwału w razie potrzeby.
• Kuncir's Algorithm 103 (1962) jest oryginalnym rekurencyjnym, dwusiecznym, adaptacyjnym integratorem. Algorytm 103 składa się z większej procedury z zagnieżdżonym podprogramem (pętla AA), która stała się rekurencyjna za pomocą goto . Zabezpiecza przed niedopełnieniem szerokości interwału (pętla BB) i przerywa, gdy tylko zostanie przekroczony określony przez użytkownika eps. Kryterium zakończenia to $, b) | <2 ^ {- n} \ epsilon$ a \ } a S ⁽²⁾ jest dokładniejszym oszacowaniem.
• Algorytm McKeemana 145 (1962) jest podobnie rekurencyjnym integratorem, który dzieli przedział na trzy zamiast na dwie części. Rekurencja jest napisana w bardziej znany sposób. Algorytm z 1962 roku, który okazał się zbyt ostrożny, używa ${\ Displaystyle | S ^ {(3)} (a, b) -S (a, b) | <3 ^ {- n} \ epsilon \,} do zakończenia, więc ulepszenie z 1963 r.$ Wykorzystuje ${\ Displaystyle {\ sqrt {3}} ^ {\, -n} \ epsilon}$ zamiast.
• Lyness (1969) jest prawie obecnym integratorem. Stworzony jako zestaw czterech modyfikacji McKeemana 1962, zastępuje trysekcję bisekcją w celu obniżenia kosztów obliczeniowych (Modyfikacje 1 + 2, zbieżne z integratorem Kuncir) i poprawia oszacowania błędów McKeemana 1962/63 do piątego rzędu (Modyfikacja 3), w sposób związany z regułą Boole'a i metodą Romberga . Modyfikacja 4, nie zaimplementowana tutaj, zawiera zapisy na błąd zaokrąglenia, który pozwala na podniesienie ε do minimum dozwolonego przez obecną precyzję i zwrócenie nowego błędu.

Linki zewnętrzne

• Moduł adaptacyjnej reguły Simpsona