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)\) приблизно скасовують. Цей безкоштовний обід, так би мовити, називається екстраполяцією Річардсона.