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

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

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

куда это интервал со средней точкой , , , и - оценки по правилу Симпсона на соответствующих интервалах, а - желаемый допуск для интервала.

Правило Симпсона является интерполяционным квадратурным правилом, которое является точным, когда подынтегральное выражение является полиномом третьей степени или ниже. С помощью Экстраполяция Ричардсона, более точная оценка Симпсона для шести значений функции сочетается с менее точной оценкой для трех значений функции путем применения поправки . Полученная таким образом оценка точна для многочленов пятой и меньшей степени.

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

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

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

  • Пусть начальный интервал будет [А, B]. Пусть исходный допуск будет ε0.
  • Для каждого подынтервала [а, б], определять D(а, б), оценка ошибки, как . Определять E = 180 ε0 / (B - А). Тогда первоначальные критерии прекращения стали бы DE.
  • Если D(а, м) ≥ D(а, б), либо достигнут уровень округления, либо ноль для ж(4) находится в интервале. Изменение толерантности ε0 к ε '0 необходимо.
    • Рекурсивные процедуры теперь должны возвращать D уровень для текущего интервала. Обычная статическая переменная E ' = 180 ε '0 / (B - А) определен и инициализирован как E.
    • (Модификация 4 i, ii) Если дальнейшая рекурсия используется на интервале:
      1. Если кажется, что округление достигнуто, измените E ' к D(а, м).[а]
      2. В противном случае отрегулируйте E ' к Максимум(E, D(а, м)).
    • Необходим некоторый контроль над настройками. Значительные увеличения и незначительные уменьшения допусков должны быть запрещены.
  • Для расчета эффективного ε '0 на всем интервале:
    • Зарегистрируйте каждый Икся на котором E ' превращается в массив (Икся, εя' ) пары. Первая запись должна быть (а, ε '0).
    • Настоящий εэфф это среднее арифметическое всех ε '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 будет превышено, заданное пользователем. Критерии прекращения действия: , куда п это текущий уровень рекурсии и S(2) это более точная оценка.[1]
  • Алгоритм МакКимана 145 (1962) представляет собой аналогичный рекурсивный интегратор, который разбивает интервал на три вместо двух частей. Рекурсия написана более привычным образом.[6] Алгоритм 1962 года, который оказался чрезмерно осторожным, использует для прекращения, поэтому улучшение 1963 года использует вместо.[3]:485,487
  • Лайнесс (1969) - почти нынешний интегратор. Созданный как набор из четырех модификаций МакКимана 1962 года, он заменяет трисекцию пополам для снижения вычислительных затрат (модификации 1 + 2, совпадающие с интегратором Кунцира) и улучшает оценки ошибок МакКимана 1962/63 годов до пятого порядка (модификация 3) в способ, связанный с Правило Буля и Метод Ромберга.[3](489) Модификация 4, которая здесь не реализована, содержит положения для ошибки округления, которая позволяет поднять ε до минимума, разрешенного текущей точностью, и вернуть новую ошибку.[3]

Примечания

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

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

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

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