Адаптивный метод Симпсонов - Adaptive Simpsons method

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

Критерий определения, когда следует прекратить деление интервала, предложенный Дж. Лайнесс,^[3] является

{ Displaystyle | S (a, m) + S (m, b) -S (a, b) | <15 epsilon ,}

куда ${ Displaystyle [а, Ь] , !}$ это интервал со средней точкой ${ Displaystyle м , !}$ , ${ Displaystyle S (а, Ь) , !}$ , ${ Displaystyle S (а, м) , !}$ , и ${ Displaystyle S (м, Ь) , !}$ - оценки по правилу Симпсона на соответствующих интервалах, а ${ Displaystyle epsilon , !}$ - желаемый допуск для интервала.

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

Численное рассмотрение

Некоторые входные данные не могут быстро сходиться в адаптивном методе Симпсона, что приводит к допуску недостаточный и создавая бесконечный цикл. Простые методы защиты от этой проблемы включают добавление ограничения глубины (как в примере C и в McKeeman), проверяя, что $ε /2 \neq ε$ в арифметике с плавающей запятой или в обоих случаях (например, Kuncir). Размер интервала также может приближаться к местному машина эпсилон, давая $а = б$ .

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

Пусть начальный интервал будет $[А, B]$ . Пусть исходный допуск будет $ε 0$ .
Для каждого подынтервала $[а, б]$ , определять $D (а, б)$ , оценка ошибки, как ${ displaystyle { frac {12} {ab}} [S ^ {(2)} (a, b) -S (a, b)] = { frac {1} {4}} [f (a) -4f (lm) + 6f (m) -4f (rm) + f (b)]}$ . Определять $E = 180 ε 0 / (B - А)$ . Тогда первоначальные критерии прекращения стали бы $D \leq E$ .
Если D(а, м) ≥ D(а, б), либо достигнут уровень округления, либо ноль для ж⁽⁴⁾ находится в интервале. Изменение толерантности ε₀ к ε '₀ необходимо.
- Рекурсивные процедуры теперь должны возвращать D уровень для текущего интервала. Обычная статическая переменная $E ' = 180 ε ' 0 / (B - А)$ определен и инициализирован как E.
- (Модификация 4 i, ii) Если дальнейшая рекурсия используется на интервале:
  1. Если кажется, что округление достигнуто, измените E ' к $D (а, м)$ .^[а]
  2. В противном случае отрегулируйте E ' к $Максимум(E, D (а, м))$ .
- Необходим некоторый контроль над настройками. Значительные увеличения и незначительные уменьшения допусков должны быть запрещены.
Для расчета эффективного ε '₀ на всем интервале:
- Зарегистрируйте каждый $Икс я$ на котором E ' превращается в массив $(Икс я, ε я')$ пары. Первая запись должна быть $(а, ε ' 0)$ .
- Настоящий ε_эфф это среднее арифметическое всех ε '₀, взвешенные по ширине интервалов.
Если нынешний E ' для интервала больше чем E, то ускорение / коррекция пятого порядка не применяется:^[b]
- Фактор «15» в критериях завершения отключен.
- Поправочный член использовать не следует.

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

Примеры реализации

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

Python

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

из __будущее__ импорт разделение # python 2 compat# "структурированная" адаптивная версия в переводе с Racketdef _quad_simpsons_mem(ж, а, фа, б, fb):    "" "Оценивает правило Симпсона, также возвращая m и f (m) для повторного использования" ""    м = (а + б) / 2    FM = ж(м)    возвращаться (м, FM, пресс(б - а) / 6 * (фа + 4 * FM + fb))def _quad_asr(ж, а, фа, б, fb, eps, весь, м, FM):    """    Эффективная рекурсивная реализация адаптивного правила Симпсона.    Значения функций в начале, середине и конце интервалов сохраняются.    """    lm, flm, оставили  = _quad_simpsons_mem(ж, а, фа, м, FM)    rm, FRM, верно = _quad_simpsons_mem(ж, м, FM, б, fb)    дельта = оставили + верно - весь    если пресс(дельта) <= 15 * eps:        возвращаться оставили + верно + дельта / 15    возвращаться _quad_asr(ж, а, фа, м, FM, eps/2, оставили , lm, flm) +\           _quad_asr(ж, м, FM, б, fb, eps/2, верно, rm, FRM)def quad_asr(ж, а, б, eps):    "" "Интегрируйте f от a до b, используя адаптивное правило Симпсона с максимальной ошибкой eps." ""    фа, fb = ж(а), ж(б)    м, FM, весь = _quad_simpsons_mem(ж, а, фа, б, fb)    возвращаться _quad_asr(ж, а, фа, б, fb, eps, весь, м, FM)из математика импорт грехРаспечатать(quad_asr(грех, 0, 1, 1e-09))

C

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

#включают  // включаем файл для fabs и sin#включают  // включаем файл для printf и perror#включают <errno.h>/ ** Адаптивное правило Симпсона, рекурсивное ядро * /плавать адаптивныйSimpsonsAux(плавать (*ж)(плавать), плавать а, плавать б, плавать eps,                          плавать весь, плавать фа, плавать fb, плавать FM, int rec) {    плавать м   = (а + б)/2,  час   = (б - а)/2;    плавать lm  = (а + м)/2,  rm  = (м + б)/2;    // серьезная числовая проблема: не сходится    если ((eps/2 == eps) || (а == lm)) { errno = ЕДОМ; возвращаться весь; }    плавать flm = ж(lm),      FRM = ж(rm);    плавать оставили  = (час/6) * (фа + 4*flm + FM);    плавать верно = (час/6) * (FM + 4*от + fb);    плавать дельта = оставили + верно - весь;    если (rec <= 0 && errno != ЕДОМ) errno = ERANGE;  // предел глубины слишком мал    // Лайнесс 1969 + экстраполяция Ричардсона; см. статью    если (rec <= 0 || фабрики(дельта) <= 15*eps)        возвращаться оставили + верно + (дельта)/15;    возвращаться адаптивныйSimpsonsAux(ж, а, м, eps/2, оставили,  фа, FM, flm, rec-1) +           адаптивныйSimpsonsAux(ж, м, б, eps/2, верно, FM, fb, от, rec-1);}/ ** Адаптивная оболочка правил Симпсона * (заполняет кешированные оценки функций) * /плавать адаптивныйSimpsons(плавать (*ж)(плавать),     // функция ptr для интеграции                       плавать а, плавать б,      // интервал [a, b]                       плавать эпсилон,         // устойчивость к ошибкам                       int maxRecDepth) {     // ограничение рекурсии    errno = 0;    плавать час = б - а;    если (час == 0) возвращаться 0;    плавать фа = ж(а), fb = ж(б), FM = ж((а + б)/2);    плавать S = (час/6)*(фа + 4*FM + fb);    возвращаться адаптивныйSimpsonsAux(ж, а, б, эпсилон, S, фа, fb, FM, maxRecDepth);}/ ** пример использования * /#включают  // для враждебного примера (функция rand)статический int callcnt = 0;статический плавать sinfc(плавать Икс) { callcnt++; возвращаться грех(Икс); } статический плавать frand48c(плавать Икс) { callcnt++; возвращаться drand48(); } int главный() {    // Пусть I будет интегралом sin (x) от 0 до 2    плавать я = адаптивныйSimpsons(sinfc, 0, 2, 1e-5, 3);    printf("интегрировать (sinf, 0, 2) =% lf п", я);   // выводим результат    перрор("adaptiveSimpsons");                   // Было ли это успешно? (глубина = 1 слишком мала)    printf("(% d оценок) п", callcnt);    callcnt = 0; srand48(0);    я = адаптивныйSimpsons(frand48c, 0, 0.25, 1e-5, 25); // случайная функция    printf("интегрировать (frand48, 0, 0,25) =% lf п", я);    перрор("adaptiveSimpsons");                        // не сходится    printf("(% d оценок) п", callcnt);    возвращаться 0;}

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

Ракетка

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

;; -----------------------------------------------------------------------------;; интерфейс, с контрактом (предоставить / контракт [адаптивный-симпсон   (-> я ((ж (-> настоящий? настоящий?)) (L настоящий?) (р  (L) (и / c настоящий? (> / c L))))       (#: эпсилон (ε настоящий?))       (р настоящий?))]);; -----------------------------------------------------------------------------;; выполнение (определять (адаптивный-симпсон ж L р #: эпсилон [ε .000000001])  (определять f @ L (ж L))  (определять f @ R (ж р))  (определение значений (M f @ M весь) (Симпсон-1call-to-f ж L f @ L р f @ R))  (asr ж L f @ L р f @ R ε весь M f @ M));; «эффективная» реализация(определять (asr ж L f @ L р f @ R ε весь M f @ M)  (определение значений (leftM  f @ leftM  оставили*)  (Симпсон-1call-to-f ж L f @ L M f @ M))  (определение значений (rightM f @ rightM верно*) (Симпсон-1call-to-f ж M f @ M р f @ R))  (определять дельта * (- (+ оставили* верно*) весь))  (cond    [(<= (пресс дельта *) (* 15 ε)) (+ оставили* верно* (/ дельта * 15))]    [еще (определять epsilon1 (/ ε 2))          (+ (asr ж L f @ L M f @ M epsilon1 оставили*  leftM  f @ leftM)              (asr ж M f @ M р f @ R epsilon1 верно* rightM f @ rightM))]));; оценить половину интервала (1 функция eval)(определять (Симпсон-1call-to-f ж L f @ L р f @ R)  (определять M (середина L р))  (определять f @ M (ж M))  (значения M f @ M (* (/ (пресс (- р L)) 6) (+ f @ L (* 4 f @ M) f @ R))))(определять (середина L р) (/ (+ L р) 2.))

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

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

Примечания

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

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

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

внешняя ссылка

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

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

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

[kuncir-1] а ^б Г.Ф. Kuncir (1962), "Алгоритм 103: интегратор правила Симпсона", Коммуникации ACM, 5 (6): 347, Дои:10.1145/367766.368179

[henriksson-2] а ^б Для более раннего нерекурсивного адаптивного интегратора, больше напоминающего Решатели ODE, видеть С. Хенрикссон (1961), "Вклад № 2: численное интегрирование Симпсона с переменной длиной шага", BIT вычислительная математика, 1: 290

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

[ornl-6] а ^б Беррилл, Марк А. "RayTrace-miniapp: src / AtomicModel / interp.hpp · de5e8229bccf60ae5c1c5bab14f861dc0326d5f9". ORNL GitLab.

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

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

[1]

[2]

[3]

[а]

[b]

[4]

[5]

[6]