3: Інтеграція

Ми хочемо побудувати числові алгоритми, які можуть виконувати певні інтеграли виду

$$I=\int_{a}^{b} f(x) d x . \nonumber$$

Обчислення цих певних інтегралів чисельно називається числовим інтегруванням, числовою квадратурою, або простіше кажучи квадратурою.

Елементарні формули

Спочатку ми розглядаємо інтеграцію від 0 до$h$, з$h$ малими, щоб служити складовими для інтеграції над більшими доменами. Ми тут визначаємо$I_{h}$ як такий інтеграл:

$$I_{h}=\int_{0}^{h} f(x) d x . \nonumber$$

Щоб виконати цей інтеграл, ми розглянемо розширення серії Тейлора$f(x)$ про значення$x=h / 2$:

$$\begin{aligned} f(x)=f(h / 2)+(x-h / 2) f^{\prime}(h / 2) &+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ +& \frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots \end{aligned} \nonumber$$

Правило середньої точки

Правило середньої точки використовує лише перший термін у розширенні серії Тейлора. Тут ми визначимо похибку в цьому наближенні. Інтеграція,

$$\begin{aligned} I_{h}=h f(h / 2)+\int_{0}^{h}((x-&h / 2) f^{\prime}(h / 2)+\frac{(x-h / 2)^{2}}{2} f^{\prime \prime}(h / 2) \\ &\left.+\frac{(x-h / 2)^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{(x-h / 2)^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d x . \end{aligned} \nonumber$$

Зміна змінних, дозволяючи$y=x-h / 2$ і$d y=d x$, і спрощення інтеграла в залежності від того, чи є integrand парним або непарним, у нас є

$$\begin{aligned} I_{h}=& h f(h / 2) \\ &+\int_{-h / 2}^{h / 2}\left(y f^{\prime}(h / 2)+\frac{y^{2}}{2} f^{\prime \prime}(h / 2)+\frac{y^{3}}{6} f^{\prime \prime \prime}(h / 2)+\frac{y^{4}}{24} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y \\ =& h f(h / 2)+\int_{0}^{h / 2}\left(y^{2} f^{\prime \prime}(h / 2)+\frac{y^{4}}{12} f^{\prime \prime \prime \prime}(h / 2)+\ldots\right) d y . \end{aligned} \nonumber$$

Інтеграли, які нам тут потрібні, є

$$\int_{0}^{\frac{h}{2}} y^{2} d y=\frac{h^{3}}{24}, \quad \int_{0}^{\frac{h}{2}} y^{4} d y=\frac{h^{5}}{160} . \nonumber$$

Тому,

$$I_{h}=h f(h / 2)+\frac{h^{3}}{24} f^{\prime \prime}(h / 2)+\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber$$

трапецієподібне правило

З розширення серії Тейлора$f(x)$ про$x=h / 2$, ми маємо

$$f(0)=f(h / 2)-\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)-\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber$$

і

$$f(h)=f(h / 2)+\frac{h}{2} f^{\prime}(h / 2)+\frac{h^{2}}{8} f^{\prime \prime}(h / 2)+\frac{h^{3}}{48} f^{\prime \prime \prime}(h / 2)+\frac{h^{4}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Додавання і множення на$h / 2$ отримуємо

$$\frac{h}{2}(f(0)+f(h))=h f(h / 2)+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Тепер підставляємо перший член праворуч, використовуючи формулу правила середньої точки:

$$\begin{aligned} \frac{h}{2}(f(0)+f(h))=\left(I_{h}-\frac{h^{3}}{24} f^{\prime \prime}(h / 2)-\frac{h^{5}}{1920} f^{\prime \prime \prime \prime}(h / 2)\right) &\, \\ &+\frac{h^{3}}{8} f^{\prime \prime}(h / 2)+\frac{h^{5}}{384} f^{\prime \prime \prime \prime}(h / 2)+\ldots, \end{aligned} \nonumber$$

і вирішуючи для$I_{h}$, знаходимо

$$I_{h}=\frac{h}{2}(f(0)+f(h))-\frac{h^{3}}{12} f^{\prime \prime}(h / 2)-\frac{h^{5}}{480} f^{\prime \prime \prime \prime}(h / 2)+\ldots \nonumber$$

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

Щоб отримати правило Сімпсона, ми об'єднаємо середню точку та трапецієподібне правило, щоб усунути термін похибки, пропорційний$h^{3}$. Помноживши (3.4) на два і додавши до (3.8), отримуємо

$$3 I_{h}=h\left(2 f(h / 2)+\frac{1}{2}(f(0)+f(h))\right)+h^{5}\left(\frac{2}{1920}-\frac{1}{480}\right) f^{\prime \prime \prime \prime}(h / 2)+\ldots, \nonumber$$

або

$$I_{h}=\frac{h}{6}(f(0)+4 f(h / 2)+f(h))-\frac{h^{5}}{2880} f^{\prime \prime \prime \prime}(h / 2)+\ldots . \nonumber$$

Зазвичай правило Сімпсона пишеться з урахуванням трьох послідовних пунктів$0, h$ і$2 h$. Підставивши$h \rightarrow 2 h$, отримуємо стандартний результат

$$I_{2 h}=\frac{h}{3}(f(0)+4 f(h)+f(2 h))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(h)+\ldots \nonumber$$

Складені правила

Тепер ми використовуємо наші елементарні формули, отримані для (3.2) для виконання інтеграла, заданого (3.1).

трапецієподібне правило

Ми припускаємо, що функція$f(x)$ відома в$n+1$ точках, позначені як$x_{0}, x_{1}, \ldots, x_{n}$, з кінцевими точками, заданими$x_{0}=a$ і$x_{n}=b$. Визначте

$$f_{i}=f\left(x_{i}\right), \quad h_{i}=x_{i+1}-x_{i} . \nonumber$$

Тоді інтеграл (3.1) може бути розкладений як

$$\begin{aligned} \int_{a}^{b} f(x) d x &=\sum_{i=0}^{n-1} \int_{x_{i}}^{x_{i+1}} f(x) d x \\ &=\sum_{i=0}^{n-1} \int_{0}^{h_{i}} f\left(x_{i}+s\right) d s, \end{aligned} \nonumber$$

де остання рівність виникає внаслідок зміни змінних$s=x-x_{i}$. Застосовуючи трапецієподібне правило до інтегралу, маємо

$$\int_{a}^{b} f(x) d x=\frac{1}{2} \sum_{i=0}^{n-1} h_{i}\left(f_{i}+f_{i+1}\right) . \nonumber$$

Якщо точки розташовані не рівномірно, скажімо, тому що дані є експериментальними значеннями, то$h_{i}$ може відрізнятися для кожного значення$i$ і (3,13) буде використовуватися безпосередньо.

Однак, якщо точки розташовані рівномірно, скажімо, тому що$f(x)$ можна обчислити, ми маємо$h_{i}=h$, незалежно від$i$. Потім ми можемо визначити

$$x_{i}=a+i h, \quad i=0,1, \ldots, n ; \nonumber$$

і оскільки кінцева точка$b$ задовольняє$b=a+n h$, ми маємо

$$h=\frac{b-a}{n} . \nonumber$$

Складене трапецієподібне правило для рівномірного простору точок потім стає

$$\begin{align} \nonumber \int_{a}^{b} f(x) d x &=\frac{h}{2} \sum_{i=0}^{n-1}\left(f_{i}+f_{i+1}\right) \\ &=\frac{h}{2}\left(f_{0}+2 f_{1}+\cdots+2 f_{n-1}+f_{n}\right) . \end{align} \nonumber$$

Перший і останній члени мають кратні одному; всі інші терміни кратні двом; і вся сума множиться на$h / 2$

Малюнок 3.1: Адаптивна квадратура Сімпсона: Рівень$1 .$

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

Тут ми розглянемо складене правило Сімпсона для рівномірних точок простору. Ми застосовуємо правило Сімпсона над інтервалами$2 h$, починаючи від$a$ і закінчуючи на$b$:

$$\begin{aligned} \int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+f_{2}\right)+\frac{h}{3}\left(f_{2}+4 f_{3}+f_{4}\right)+\ldots & \\ &+\frac{h}{3}\left(f_{n-2}+4 f_{n-1}+f_{n}\right) . \end{aligned} \nonumber$$

Зверніть увагу, що$n$ повинна бути навіть, щоб ця схема працювала. Поєднуючи терміни, ми маємо

$$\int_{a}^{b} f(x) d x=\frac{h}{3}\left(f_{0}+4 f_{1}+2 f_{2}+4 f_{3}+2 f_{4}+\cdots+4 f_{n-1}+f_{n}\right) . \nonumber$$

Перший і останній терміни кратні одному; парні індексовані терміни кратні 2; непарні індексовані терміни кратні 4; і вся сума множиться на$h / 3$.

Адаптивна інтеграція

Корисна функція MATLAB quad.m виконує числове інтегрування з використанням адаптивної квадратури Сімпсона. Ідея полягає в тому, щоб дозволити обчисленню самому визначитися з розміром сітки, необхідним для досягнення певного рівня точності. Більше того, розмір сітки не повинен бути однаковим у всьому регіоні інтеграції.

Ми починаємо адаптивну інтеграцію на тому, що називається Level 1. Рівномірно розташовані точки, в яких повинна бути$f(x)$ оцінена функція, показані на рис.3.1. Відстань між точками$a$ і$b$ приймається бути$2 h$, щоб

$$h=\frac{b-a}{2} . \nonumber$$

Інтеграція за допомогою правила Сімпсона (3.11) з розміром сітки$h$ дає інтеграл$I$,

$$I=\frac{h}{3}(f(a)+4 f(c)+f(b))-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \nonumber$$

де$\xi$ деяка цінність задовольняє$a \leq \xi \leq b$. Інтеграція за допомогою правила Сімпсона двічі з$h / 2$ прибутковістю розміру сітки

$$I=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b))-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{l}\right)-\frac{(h / 2)^{5}}{90} f^{\prime \prime \prime \prime}\left(\xi_{r}\right), \nonumber$$

з$\xi_{l}$ і$\xi_{r}$ деякими цінностями задовольняють$a \leq \xi_{l} \leq c$ і$c \leq \xi_{r} \leq b$.

Тепер ми визначимо два наближення до інтеграла

$$\begin{aligned} &S_{1}=\frac{h}{3}(f(a)+4 f(c)+f(b)), \\ &S_{2}=\frac{h}{6}(f(a)+4 f(d)+2 f(c)+4 f(e)+f(b)), \end{aligned} \nonumber$$

і дві пов'язані помилки

$$\begin{aligned} &E_{1}=-\frac{h^{5}}{90} f^{\prime \prime \prime \prime}(\xi), \\ &E_{2}=-\frac{h^{5}}{2^{5} \cdot 90}\left(f^{\prime \prime \prime \prime}\left(\xi_{l}\right)+f^{\prime \prime \prime \prime}\left(\xi_{r}\right)\right) . \end{aligned} \nonumber$$

Тепер ми запитуємо, чи достатньо точне значення$S_{2}$ для інтеграла, або ми повинні додатково вдосконалити обчислення і перейти до рівня 2? Щоб відповісти на це питання, ми зробимо спрощене наближення, що всі похідні четвертого порядку$f(x)$ в термах похибки рівні; тобто

$$f^{\prime \prime \prime \prime}(\xi)=f^{\prime \prime \prime \prime}\left(\xi_{l}\right)=f^{\prime \prime \prime \prime}\left(\xi_{r}\right)=C . \nonumber$$

Тоді

$$\begin{aligned} &E_{1}=-\frac{h^{5}}{90} C, \\ &E_{2}=-\frac{h^{5}}{2^{4} \cdot 90} C=\frac{1}{16} E_{1} . \end{aligned} \nonumber$$

Тепер, оскільки інтеграл дорівнює наближенню плюс пов'язана з ним похибка,

$$S_{1}+E_{1}=S_{2}+E_{2} \text {, } \nonumber$$

і з тих пір

$$E_{1}=16 E_{2} \text {, } \nonumber$$

ми можемо вивести оцінку терміна помилки$E_{2}$:

$$E_{2}=\frac{1}{15}\left(S_{2}-S_{1}\right) \text {. } \nonumber$$

Тому, враховуючи деяку специфічну величину допуску тол, якщо

$$\left|\frac{1}{15}\left(S_{2}-S_{1}\right)\right|<\text { tol, } \nonumber$$

тоді ми можемо прийняти$S_{2}$ як$I$. Якщо оцінка похибки більше за величиною, ніж tol, то переходимо до рівня 2. Обчислення на рівні 2 додатково ділить інтервал інтеграції від$a$ до на два інтервали інтеграції$a$ до$c$ і$c$ до$b$, і$b$ продовжується з вищевказаною процедурою. самостійно на обох половинках. Інтеграція може бути зупинена на будь-якій половині за умови, що допуск менше tol/2 (оскільки сума обох помилок повинна бути меншою за tol). В іншому випадку будь-яка половина може перейти до рівня 3, і так далі.

Як примітка сторони, два значення$I$ наведені вище (для інтеграції з розміром кроку$h$ і$h / 2)$ можуть бути об'єднані, щоб дати більш точне значення для I, заданого

$$I=\frac{16 S_{2}-S_{1}}{15}+\mathrm{O}\left(h^{7}\right), \nonumber$$

де терміни помилки$\mathrm{O}\left(h^{5}\right)$ приблизно скасовують. Цей безкоштовний обід, так би мовити, називається екстраполяцією Річардсона.