Адаптивный метод Симпсона

Адаптивный метод Симпсона , также называемый адаптивным правилом Симпсона , — это метод численного интегрирования, предложенный Г. Ф. Кунциром в 1962 году. ^{[ 1 ]} Вероятно, это первый рекурсивный адаптивный алгоритм численного интегрирования, появившийся в печати. ^{[ 2 ]} более современные адаптивные методы, основанные на квадратуре Гаусса – Кронрода и квадратуре Кленшоу – Кертиса хотя сейчас обычно предпочитаются . Адаптивный метод Симпсона использует оценку ошибки, которую мы получаем при вычислении определенного интеграла с использованием правила Симпсона . Если ошибка превышает заданный пользователем допуск, алгоритм требует разделить интервал интегрирования на два и рекурсивно применить адаптивный метод Симпсона к каждому подинтервалу. Этот метод обычно намного более эффективен, чем составное правило Симпсона , поскольку он использует меньше вычислений функции в тех местах, где функция хорошо аппроксимируется кубической функцией .

Правило Симпсона — это интерполяционное правило квадратур, которое является точным, когда подынтегральная функция является многочленом третьей степени или ниже. Используя экстраполяцию Ричардсона , можно получить более точную оценку Симпсона. $S(a,m)+S(m,b)$ для шести значений функции сочетается с менее точной оценкой $S(a,b)$ для трех значений функции, применяя поправку $[S(a,m)+S(m,b)-S(a,b)]/15$ . Таким образом, полученная оценка точна для полиномов пятой и меньшей степени.

Математическая процедура

Определение терминов

Критерий определения момента прекращения деления интервала, предложенный Дж. Н. Лайнсом, ^{[ 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 и в Маккимане), проверку того, что $ε /2 \neq ε$ в арифметике с плавающей запятой, или и то, и другое (как в Кунцире). Размер интервала также может приближаться к эпсилону локальной машины , что дает $a = b$ .

Статья Лайнса 1969 года включает «Модификацию 4», которая решает эту проблему более конкретным образом: ^{[ 3 ]}^: 490–2

Пусть начальный интервал будет $[A, B]$ . Пусть исходный допуск равен $ε 0$ .
Для каждого подинтервала $[a, b]$ определите $D (a, b)$ , оценку ошибки, как ${\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)]}$ . Определить $E знак равно 180 ε 0 / (B - A)$ . исходными критериями завершения станут $D \leq E.$ Тогда
Если D ( a , m ) ≥ D ( a , b ) , либо достигнут уровень округления, либо нуль для f ⁽⁴⁾ находится в интервале. изменить допуск ε ₀ на ε′ ₀ Необходимо .
- Рекурсивным процедурам теперь необходимо возвращать уровень D для текущего интервала. Подпрограммно-статическая переменная $E' = 180 ε' 0 / (B - A)$ определяется и инициализируется значением E .
- (Модификация 4 i, ii) Если на интервале используется дальнейшая рекурсия:
  1. Если кажется, что округление достигнуто, измените E' на $D (a, m)$ . ^{[ а ]}
  2. В противном случае настройте E' на $max(E, D (a, m))$ .
- Необходим некоторый контроль над корректировками. Значительное увеличение и незначительное уменьшение допусков следует запрещать.
Чтобы вычислить эффективное ε′ ₀ на всем интервале:
- Зарегистрируйте каждый $x i,$ в котором E' изменяется в массив $пар (x i, ε i')$ . Первая запись должна быть $(a, ε' 0)$ .
- Фактическое ε _eff представляет собой среднее арифметическое всех ε' ₀ , взвешенное по ширине интервалов.
Если текущее E' для интервала выше, чем E , то ускорение/коррекция пятого порядка не будет применяться: ^{[ б ]}
- Коэффициент «15» в критериях завершения отключен.
- Поправочный термин использовать не следует.

Маневр повышения эпсилона позволяет использовать программу в режиме «максимальных усилий»: при нулевом начальном допуске программа попытается получить наиболее точный ответ и вернуть фактический уровень ошибки. ^{[ 3 ]}^: 492

Примеры реализации кода

Общий метод реализации, показанный ниже, — это передача $f(a), f(b), f(m)$ вместе с интервалом $[a, b]$ . Эти значения, используемые для оценки $S (a, b)$ на родительском уровне, снова будут использоваться для подинтервалов. Это сокращает стоимость каждого рекурсивного вызова с 6 до 2 вычислений входной функции. Размер используемого стекового пространства остается линейным в зависимости от уровня рекурсий.

Питон

Вот реализация адаптивного метода Симпсона на Python .

from __future__ import division # python 2 compat
# "structured" adaptive version, translated from Racket
def _quad_simpsons_mem(f, a, fa, b, fb):
    """Evaluates the Simpson's Rule, also returning m and f(m) to reuse"""
    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):
    """
    Efficient recursive implementation of adaptive Simpson's rule.
    Function values at the start, middle, end of the intervals are retained.
    """
    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, right, rm, frm)

def quad_asr(f, a, b, eps):
    """Integrate f from a to b using Adaptive Simpson's Rule with max error of 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, whole, m, fm)

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

С

Вот реализация адаптивного метода Симпсона в C99, которая позволяет избежать избыточных оценок f и квадратурных вычислений. Он включает в себя все три «простых» способа защиты от числовых проблем.

#include <math.h>  // include file for fabs and sin
#include <stdio.h> // include file for printf and perror
#include <errno.h>

/** Adaptive Simpson's Rule, Recursive Core */
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;
    // serious numerical trouble: it won't converge
    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);
    float right = (h/6) * (fm + 4*frm + fb);
    float delta = left + right - whole;

    if (rec <= 0 && errno != EDOM) errno = ERANGE;  // depth limit too shallow
    // Lyness 1969 + Richardson extrapolation; see article
    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);
}

/** Adaptive Simpson's Rule Wrapper
 *  (fills in cached function evaluations) */
float adaptiveSimpsons(float (*f)(float),     // function ptr to integrate
                       float a, float b,      // interval [a,b]
                       float epsilon,         // error tolerance
                       int maxRecDepth) {     // recursion cap
    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);
}

/** usage example */
#include <stdlib.h> // for the hostile example (rand function)
static int callcnt = 0;
static float sinfc(float x) { callcnt++; return sinf(x); } 
static float frand48c(float x) { callcnt++; return drand48(); } 
int main() {
    // Let I be the integral of sin(x) from 0 to 2
    float I = adaptiveSimpsons(sinfc, 0, 2, 1e-5, 3);
    printf("integrate(sinf, 0, 2) = %lf\n", I);   // print the result
    perror("adaptiveSimpsons");                   // Was it successful? (depth=1 is too shallow)
    printf("(%d evaluations)\n", callcnt);

    callcnt = 0; srand48(0);
    I = adaptiveSimpsons(frand48c, 0, 0.25, 1e-5, 25); // a random function
    printf("integrate(frand48, 0, 0.25) = %lf\n", I);
    perror("adaptiveSimpsons");                        // won't converge
    printf("(%d evaluations)\n", callcnt);
    return 0;
}

Эта реализация была включена в программу трассировки лучей C++, предназначенную для моделирования рентгеновского лазера в Национальной лаборатории Ок-Ридж . ^{[ 4 ]} среди других проектов. Версия ORNL дополнена счетчиком вызовов, шаблонами для различных типов данных и оболочками для интеграции по нескольким измерениям. ^{[ 4 ]}

Ракетка

Вот реализация адаптивного метода Симпсона в Racket с поведенческим программным контрактом. Экспортированная функция вычисляет неопределенный интеграл для некоторой заданной функции f . ^{[ 5 ]}

;; -----------------------------------------------------------------------------
;; interface, with contract 
(provide/contract
 [adaptive-simpson 
  (->i ((f (-> real? real?)) (L real?) (R  (L) (and/c real? (>/c L))))
       (#:epsilon (ε real?))
       (r real?))])

;; -----------------------------------------------------------------------------
;; implementation 
(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))

;; the "efficient" implementation
(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))]))

;; evaluate half an interval (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))))

(define (mid L R) (/ (+ L R) 2.))

Связанные алгоритмы

Хенрикссон (1961) представляет собой нерекурсивный вариант правила Симпсона. Он «адаптируется», интегрируя слева направо и регулируя ширину интервала по мере необходимости. ^{[ 2 ]}
Алгоритм Кунцира 103 (1962) является оригинальным рекурсивным, делящим пополам, адаптивным интегратором. Алгоритм 103 состоит из более крупной процедуры с вложенной подпрограммой (цикл AA), выполненной рекурсивной за счет использования оператора goto . Он защищает от потери ширины интервала (цикл BB) и прерывается, как только будет превышено заданное пользователем значение eps. Критерием прекращения является $|S^{(2)}(a,b)-S(a,b)|<2^{-n}\epsilon \,$ , где $n$ — текущий уровень рекурсии, а $S (2)$ это более точная оценка. ^{[ 1 ]}
Алгоритм МакКимана 145 (1962) представляет собой аналогичный рекурсивный интегратор, который делит интервал на три, а не на две части. Рекурсия написана более привычным образом. ^{[ 6 ]} Алгоритм 1962 года, оказавшийся слишком осторожным, использует $|S^{(3)}(a,b)-S(a,b)|<3^{-n}\epsilon \,$ для завершения, поэтому в усовершенствовании 1963 года используется ${\sqrt {3}}^{\,-n}\epsilon$ вместо. ^{[ 3 ]}^{: 485, 487}
Лайнесс (1969) — почти нынешний интегратор. Созданный как набор из четырех модификаций МакКимана 1962 года, он заменяет трисекцию на биссекцию для снижения вычислительных затрат (Модификации 1+2, совпадающие с интегратором Кунцира) и улучшает оценки ошибок МакКимана 1962/63 года до пятого порядка (Модификация 3), в способ, связанный с правилом Буля и методом Ромберга . ^{[ 3 ]}^: 489 Модификация 4, не реализованная здесь, содержит положения об ошибке округления, которые позволяют повысить ε до минимума, допустимого текущей точностью, и вернуть новую ошибку. ^{[ 3 ]}

Примечания

^ В оригинальном 4i упоминается только повышение E'. Однако в более позднем тексте упоминается, что его можно снижать большими шагами.
^ Это, вероятно, также относится к уменьшению ширины допуска/интервала в упрощенном случае.

Библиография

^ Перейти обратно: ^а ^б Г. Ф. Кунчир (1962), «Алгоритм 103: интегратор правил Симпсона», Communications of ACM , 5 (6): 347, doi : 10.1145/367766.368179
^ Перейти обратно: ^а ^б Информацию о более раннем нерекурсивном адаптивном интеграторе, больше напоминающем решатели ОДУ , см. С. Хенрикссон (1961), «Вклад № 2: Численное интегрирование Симпсона с переменной длиной шага», BIT Numerical Mathematics , 1 : 290
^ Перейти обратно: ^а ^б ^с ^д ^и ^ж Дж. Н. Лайнесс (1969), «Заметки об адаптивной квадратурной программе Симпсона», Журнал ACM , 16 (3): 483–495, doi : 10.1145/321526.321537
^ Перейти обратно: ^а ^б Беррилл, Марк А. «RayTrace-miniapp: src/AtomicModel/interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9» . ОРНЛ GitLab .
^ Феллизен, Маттиас. «[рэкет] адаптивная интеграция Симпсона» . Список рассылки рэкета «Пользователи» . Проверено 26 сентября 2018 г.
^ Маккиман, Уильям Маршалл (1 декабря 1962 г.). «Алгоритм 145: Адаптивное численное интегрирование по правилу Симпсона» . Коммуникации АКМ . 5 (12): 604. дои : 10.1145/355580.369102 .

Внешние ссылки

Модуль адаптивного правила Симпсона

[4] В оригинальном 4i упоминается только повышение E'. Однако в более позднем тексте упоминается, что его можно снижать большими шагами.

[5] Это, вероятно, также относится к уменьшению ширины допуска/интервала в упрощенном случае.

[kuncir-1] Перейти обратно: ^а ^б Г. Ф. Кунчир (1962), «Алгоритм 103: интегратор правил Симпсона», Communications of ACM , 5 (6): 347, doi : 10.1145/367766.368179

[henriksson-2] Перейти обратно: ^а ^б Информацию о более раннем нерекурсивном адаптивном интеграторе, больше напоминающем решатели ОДУ , см. С. Хенрикссон (1961), «Вклад № 2: Численное интегрирование Симпсона с переменной длиной шага», BIT Numerical Mathematics , 1 : 290

[lyness-3] Перейти обратно: ^а ^б ^с ^д ^и ^ж Дж. Н. Лайнесс (1969), «Заметки об адаптивной квадратурной программе Симпсона», Журнал ACM , 16 (3): 483–495, doi : 10.1145/321526.321537

[ornl-6] Перейти обратно: ^а ^б Беррилл, Марк А. «RayTrace-miniapp: src/AtomicModel/interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9» . ОРНЛ GitLab .

[7] Феллизен, Маттиас. «[рэкет] адаптивная интеграция Симпсона» . Список рассылки рэкета «Пользователи» . Проверено 26 сентября 2018 г.

[8] Маккиман, Уильям Маршалл (1 декабря 1962 г.). «Алгоритм 145: Адаптивное численное интегрирование по правилу Симпсона» . Коммуникации АКМ . 5 (12): 604. дои : 10.1145/355580.369102 .

[ 1 ]

[ 2 ]

[ 3 ]

[ а ]

[ б ]

[ 4 ]

[ 5 ]

[ 6 ]

v т и Численное интегрирование
Newton–Cotes formulas	Trapezoidal rule Simpson's rule Simpson's 3/8 rule Adaptive Simpson's method Boole's rule Romberg's method
Gaussian quadrature	Gauss–Hermite quadrature Gauss–Jacobi quadrature Gauss–Kronrod quadrature formula Gauss–Laguerre quadrature Gauss–Legendre quadrature Chebyshev–Gauss quadrature
Other	Barnes–Hut simulation Bayesian quadrature Clenshaw–Curtis quadrature Filon quadrature Lebedev quadrature Tanh-sinh quadrature