9.3: Правило Сімпсона

Наш аналіз правила середньої точки та правила трапеції показав, що обидва методи мають$O(1/N^2)$ числову похибку. В обох випадках похибка може бути простежена до одного і того ж джерела: той факт, що інтегральна оцінка по кожному сегменту відрізняється від результату ряду Тейлора в другому порядку - термін, пропорційний$f''(x_n)$. Це пропонує спосіб покращити результат числового інтегрування: ми могли б взяти середньозважене середнє значення правила середньої точки та правила трапеції, таким чином, що числові помилки другого порядку з двох схем скасовують один одного! Це метод числового інтегрування, відомий як правило Сімпсона.

Якщо бути точним, давайте знову розглянемо пару сусідніх відрізків, які лежать між рівнорозташованими точками дискретизації$\left\{x_{n-1}, x_n, x_{n+1}\right\}$. Як виведено вище, інтеграл над цими сегментами може бути розширений Тейлором як

$$\mathcal{I}_n = \; 2 f(x_n) \Delta x \;+\; \frac{f''(x_n)}{3} \Delta x^3 + O(\Delta x^5) \;+\; \cdots$$

Для порівняння, оцінювачі правила середньої точки та правила трапеції для інтеграла є

$$\mathcal{I}_n^{\mathrm{mp}} = 2f(x_n) \Delta x$$

$$\mathcal{I}_n^{\mathrm{trapz}} = 2 f(x_n) \Delta x + \frac{f''(x_n)}{2} \Delta x^3 + O(\Delta x^5).$$

Отже, ми могли б взяти наступне середньозважене:

$$\mathcal{I}_n^{\mathrm{simp}} = \frac{1}{3}\mathcal{I}_n^{\mathrm{mp}} + \frac{2}{3}\mathcal{I}_n^{\mathrm{trapz}} = 2 f(x_n) \Delta x + \frac{f''(x_n)}{3} \Delta x^3 + O(\Delta x^5).$$

Таке середньозважене буде відповідати результату серії Тейлора аж до$O(\Delta x^4)$! (Ви можете перевірити самі, чи відрізняються$O(\Delta x^5)$ умови.) Підводячи підсумок, правило Сімпсона для цього набору з трьох пунктів можна записати як

$$\begin{align} \mathcal{I}_n^{\mathrm{simp}} &= \frac{1}{3} \Big[2f(x_n) \Delta x\Big] + \frac{2}{3} \, \Delta x \Big[ \frac{f(x_{n-1})}{2} + f(x_n) + \frac{f(x_{n+1})}{2}\Big]\\ &= \frac{\Delta x}{3} \Big[f(x_{n-1}) + 4f(x_n) + f(x_{n+1})\Big]. \end{align}$$

Загальна числова похибка, над набором$O(N)$ сегментів, дорівнює$O(1/N^4)$. Це вдосконалення двох повноважень$1/N$ над правилами трапцезіуму та середньої точки! Що ще краще, це те, що він включає в себе точно таку ж кількість арифметичних операцій, що і правило трапеції. Це так близько до безкоштовного обіду, як ви можете отримати в обчислювальній науці.

9.3.1 Реалізація правила Сімпсона на Python

У Scipy правило Сімпсона реалізовано функцією scipy.integrate.simps, яка визначається у підмодулі scipy.integrate. Подібно до функції trapz, це може бути викликано або simps (y, x) або simps (y, dx=s) для оцінки інтеграла$\int y \;dx$, використовуючи елементи x як точки дискретизації, при цьому y вказуючи набір значень для integrand.

Оскільки правило Сімпсона вимагає поділу сегментів на пари, якщо ви вкажете парну кількість точок дискретизації в x (тобто непарну кількість сегментів), функція вирішує це, виконуючи оцінку правила трапеції на першому та останньому сегментах. Зазвичай помилка незначна, тому не турбуйтеся про цю деталь

Ось приклад симпсів в дії:

>>> from scipy import *
>>> from scipy.integrate import simps
>>> x = linspace(0,10,25)
>>> y = exp(-x)
>>> t = simps(y,x)
>>> print(t)
1.00011864276

Для такої ж кількості точок дискретизації дає правило трапеції$1.01438$; точний результат -$0.9999546\dots$ Ясно, правило Сімпсона більш точне.