Формула Симпсона

Формула Симпсона (также Ньютона-Симпсона^[1]) относится к приёмам численного интегрирования. Получила название в честь британского математика Томаса Симпсона (1710—1761).

Суть метода заключается в приближении подынтегральной функции на отрезке [math]\displaystyle{ [a,b] }[/math] интерполяционным многочленом второй степени [math]\displaystyle{ p_2(x) }[/math], то есть приближение графика функции на отрезке параболой. Метод Симпсона имеет порядок погрешности 4 и алгебраический порядок точности 3.

Формула

Формулой Симпсона называется интеграл от интерполяционного многочлена второй степени на отрезке [math]\displaystyle{ [a,b] }[/math]:

[math]\displaystyle{ {\int\limits_a^b f(x) dx} \approx {\int\limits_{a}^{b} {p_2(x)} dx} = \frac{b-a}{6}{ \left( f(a) + 4 f\left(\frac{a+b}{2}\right) + f(b) \right)}, }[/math]

где [math]\displaystyle{ f(a) }[/math], [math]\displaystyle{ f((a+b)/2) }[/math] и [math]\displaystyle{ f(b) }[/math] — значения функции в соответствующих точках (на концах отрезка и в его середине).

Погрешность

При условии, что у функции [math]\displaystyle{ f(x) }[/math] на отрезке [math]\displaystyle{ [a,b] }[/math] существует четвёртая производная, погрешность [math]\displaystyle{ E(f) }[/math], согласно найденной Джузеппе Пеано формуле, равна:

[math]\displaystyle{ E(f) = - \frac{(b-a)^5}{2880}{{f^{(4)}(\zeta)}}, \ \ \ \zeta \in [a,b]. }[/math]

В связи с тем, что значение [math]\displaystyle{ \zeta }[/math] зачастую неизвестно, для оценки погрешности используется следующее неравенство:

[math]\displaystyle{ \left| E(f) \right| \leqslant \frac{(b-a)^5}{2880} \max\limits_{x\in[a,b]} {\left| f^{(4)}(x) \right|}. }[/math]

Представление в виде метода Рунге-Кутты

Формулу Симпсона можно представить в виде таблицы метода Рунге-Кутты следующим образом:

[math]\displaystyle{ \begin{array}{c|ccc} 0 &&&\\ \frac{1}{2} & \frac{1}{2} &&\\ 1 & -1 & 2&\\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} }[/math]

Составная формула (формула Котеса)

Для более точного вычисления интеграла интервал [math]\displaystyle{ [a,b] }[/math] разбивают на [math]\displaystyle{ N=2n }[/math] элементарных отрезков одинаковой длины и применяют формулу Симпсона на составных отрезках. Каждый составной отрезок состоит из соседней пары элементарных отрезков. Значение исходного интеграла является суммой результатов интегрирования на составных отрезках:

[math]\displaystyle{ \int_a^b f(x) \, dx\approx \frac{h}{3}\left[f(x_0)+2\sum_{j=1}^{N/2-1}f(x_{2j})+ 4\sum_{j=1}^{N/2}f(x_{2j-1})+f(x_N) \right], }[/math]

где [math]\displaystyle{ h = \frac{b-a}{N} }[/math] — величина шага, а [math]\displaystyle{ x_j=a+j h }[/math] — чередующиеся границы и середины составных отрезков, на которых применяется формула Симпсона. Один подобный составной отрезок [math]\displaystyle{ [x_{j-1}, x_{j+1}] }[/math]состоит из двух элементарных отрезков [math]\displaystyle{ [x_{j-1}, x_{j}], [x_{j}, x_{j+1}] }[/math]. Таким образом, если проводить параллели с простой формулой Симпсона, то в данном случае серединой отрезка, на котором применяется формула Симпсона, становится [math]\displaystyle{ x_j }[/math].

Обычно для равномерной сетки данную формулу записывают в других обозначениях (отрезок [math]\displaystyle{ [a,b] }[/math] разбит на [math]\displaystyle{ N }[/math] отрезков) в виде

[math]\displaystyle{ \int_a^b f(x) \, dx\approx \frac{h}{3}\bigg[f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+2f(x_4)+\cdots+4f(x_{N-1})+f(x_N)\bigg]. }[/math]

Также формулу можно записать используя только известные значения функции, то есть значения в узлах:

[math]\displaystyle{ {\int\limits_a^b f(x) dx} \approx \frac{h}{3} \cdot \sum_{k=1,2}^{N-1} \left[ f(x_{k-1})+4f(x_k)+f(x_{k+1}) \right] }[/math]

где [math]\displaystyle{ k=1,2 }[/math] означает что индекс меняется от единицы с шагом, равным двум.

Общая погрешность [math]\displaystyle{ E(f) }[/math] при интегрировании по отрезку [math]\displaystyle{ [a,b] }[/math] с шагом [math]\displaystyle{ x_i - x_{i-1} = h }[/math] (при этом, в частности, [math]\displaystyle{ x_0 = a }[/math], [math]\displaystyle{ x_N = b }[/math]) определяется по формуле^[2]:

[math]\displaystyle{ \left| E(f) \right| \leqslant \frac{(b-a)}{2880}h^4 \max\limits_{x\in[a,b]} |f^{(4)} (x)| }[/math].

При невозможности оценить погрешность с помощью максимума четвёртой производной (например, на заданном отрезке она не существует, либо стремится к бесконечности), можно использовать более грубую оценку:

[math]\displaystyle{ \left| E(f) \right| \leqslant \frac{(b-a)}{288}h^3 \max\limits_{x\in[a,b]} |f^{(3)} (x)| }[/math].

Проверка составной формулы Симпсона в случае интегрирования узких пиков

Составная формула Симпсона не проходит проверку на величину погрешности в случае узких (малое число точек на пик) пикоподобных функций, оказываясь значительно менее эффективной^[3], чем правило трапеций. Именно, для достижения той же погрешности, что и в случае правила трапеций, составному правилу Симпсона требуется в 1.8 раз больше точек. Интеграл по составному правилу Симпсона может быть разложен на суперпозицию двух интегралов: 2/3 интеграла по правилу трапеций с шагом h, и 1/3 интеграла по правилу центральных прямоугольников с шагом 2h, и погрешность составного правила Симпсона соответствует второму слагаемому. Можно построить удовлетворительную модификацию правила Симпсона путём усреднения схем этого правила, полученных со сдвигом рамки суммирования на одну точку, при этом получаются следующие правила^[3]:[math]\displaystyle{ \int_a^b f(x) \, dx\approx \tfrac{h}{24} \left[-f(x_{-1})+12f(x_0)+25f(x_1)+24 \sum_{i=2}^{n-2} f(x_i)+25f(x_{n-1})+12f(x_n)-f(x_{n+1})\right] }[/math]в котором используются значения, выходящие за границу интервала интегрирования, или[math]\displaystyle{ \int_a^b f(x) \, dx\approx \tfrac{h}{24} \left[9f(x_0)+28f(x_1)+23f(x_2)+24 \sum_{i=3}^{n-3} f(x_i)+23f(x_{n-2})+28f(x_{n-1})+9f(x_n)\right] }[/math]в котором значения, находящиеся за границей интервала интегрирования, не используются. Приложение второго из правил к участку из трёх точек порождает правило Симпсона 1/3, к участку из 4 точек - 3/8.

В этих правилах веса точек внутри интервала интегрирования равны единице, отличия наблюдаются только на концах участка. Эти правила могут быть ассоциированы с формулой Эйлера-Маклорена, при условии учета первой производной и названы правилами Эйлера-Маклорена первого порядка^[3]. Разница между правилами состоит в способе вычисления первой производной на краях интервала интегрирования. Разница первых производных на краях участка интегрирования учитывает вклад второй производной в интеграл функции. Формула Эйлера-Маклорена аналогично приведённым выше правилам первого порядка может быть использована для конструирования правил интегрирования третьего, пятого и более высоких порядков.

См. также

Примечания

↑ Формула Ньютона-Симпсона (неопр.) (недоступная ссылка). Дата обращения: 14 августа 2009. Архивировано 22 мая 2010 года.
↑ Численные методы / Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. — 4-е изд. — М.: БИНОМ, Лаборатория знаний, 2006. — С. 122. — 636 с. — ISBN 5-94774-396-5.
↑ ^3,0 ^3,1 ^3,2 Comparison of integration rules in the case of very narrow chromatographic peaks (англ.) // Chemometrics and Intelligent Laboratory Systems. — 2018-08-15. — Vol. 179. — P. 22–30. — ISSN 0169-7439. — doi:10.1016/j.chemolab.2018.06.001.

Литература

Костомаров Д. П., Фаворский А. П. Вводные лекции по численным методам. М.: Логос, 2004. 184 с. ISBN 5-94010-286-7
Петров И. Б., Лобанов А. И. Лекции по вычислительной математике. М.: Интуит, Бином, 2006. 523 с. ISBN 5-94774-542-9

[1] Формула Ньютона-Симпсона (неопр.) (недоступная ссылка). Дата обращения: 14 августа 2009. Архивировано 22 мая 2010 года.

[2] Численные методы / Н. С. Бахвалов, Н. П. Жидков, Г. М. Кобельков. — 4-е изд. — М.: БИНОМ, Лаборатория знаний, 2006. — С. 122. — 636 с. — ISBN 5-94774-396-5.

[:0-3] 3,0 ^3,1 ^3,2 Comparison of integration rules in the case of very narrow chromatographic peaks (англ.) // Chemometrics and Intelligent Laboratory Systems. — 2018-08-15. — Vol. 179. — P. 22–30. — ISSN 0169-7439. — doi:10.1016/j.chemolab.2018.06.001.

[1]

[2]

[3]