Подвійні та потрійні інтеграли - Практичне застосування чисельних методів мовою Python
Основні поняття. Чисельні методи лінійної алгебри та аналізу
Зміст
Швидкий пошук
Метод прямокутників для подвійного інтеграла
Розглянемо питання чисельного знаходження подвійного інтеграла по прямокутній області \([a,b] \times [c, d]\) :
Розглянемо два способи виведення аналога формули прямокутників для чисельного розв'язання інтеграла.
Висновок за допомогою одновимірних інтегралів¶
Так як ми вже знаємо формули прямокутників для інтегралів від функцій однієї змінної, то зручно уявити подвійний інтеграл через два інтеграли, кожен з яких обчислюватиметься від функції однієї змінної і може бути чисельно знайдений за допомогою вже відомої одномірної формули прямокутників. З цією метою введемо допоміжну функцію \(g(x)\):
Кожен із інтегралів
можна обчислити з допомогою формул чисельного інтегрування одномірних інтегралів. Скористаємося формулою прямокутників ( 3.1.5 ) і почнемо з \(g(x) = \int_c^d f(x, y) dy\) . На відрізку \([c, d]\) введемо рівномірну сітку з кроком \(h_y\) :
Тоді для аналізованого інтеграла складова формула прямокутників матиме вигляд:
Для чисельного знаходження інтеграла по змінній \(x\) на відрізку \([a, b]\) введемо рівномірну сітку з кроком \(h_y\) :
на якій для інтеграла \(\int_a^b g(x) dx\) отримаємо формулу прямокутників:
Об'єднуючи формули в кожному напрямку, отримаємо складову формулу прямокутників для подвійного інтеграла:
Прямий висновок
Формула (1) може бути отримана безпосередньо з використанням основної ідеї побудови одновимірної формули прямокутників. Розіб'ємо прямокутник \([a, b] \times [c, d]\) на \(N_x \times N_y\) осередків(Прямокутників). Ідея методу прямокутників полягає наближенні функції \(f\) шматково-постійною функцією: з постійними значеннями на осередках рівними значенню \(f\) в центрі осередку. Осередок \((i, j)\) займає область
центром якої є точка ((x_i, y_j) \in \omega_h = \omega_ \times \omega_\) . Таким чином, інтеграл по комірці дорівнює \(h_xh_y f(x_i, y_j)\) , а подвійний інтеграл по всій області є сума по всіх комірках, що дає формулу (1).
Програмна реалізація подвійної суми¶
Формула (1) містить подвійну суму, яку реалізують як подвійний цикл for . Функція, що реалізує формулу (1) мовою Python, може виглядати так:
Зберігши цю функцію у файлі-модулі rectangular_double.py , ми можемо обчислити інтеграл, наприклад, \(\int_2^3\int_0^2 (2x + y) dxdy\) в інтерактивній оболонці і показати, що функція повертає правильне значення:
Повторне використання коду для одновимірного інтеграла¶
Цілком природно реалізувати метод прямокутників для чисельного знаходження подвійного інтеграла оскільки ми зробили функції rectangular_double1 за умови, що ми знаємо формулу (1). Проте, постає питання: чи можемо ми (як це робили при виведенні формули через формулу для одномірного інтеграла) повторно використати реалізацію для одномірного інтеграла при обчисленні подвійного? Тобто, чи можемо ми використовувати функцію rectangular із розділу Реалізація “двічі”? Відповідь: "Так, можемо". Відповідна функція буде дуже короткою:
Важлива перевага такого підходу полягає в тому, що ми використовуємо вже протестовану функцію для одномірного методу прямокутників і використовуємо одномірний інтеграл для обчислення двовимірного точно як і робили це математично.
Перевірка за допомогоютестових функцій¶
Як можна випробувати, що функція для подвійного інтеграла працює правильно? Кращий юніт тест — знайти завдання, де похибка дорівнювала нулю, тому що в цьому випадку ми точно знатимемо чисельну відповідь. Одновимірна формула прямокутників точна на лінійних функціях незалежно від кількості часткових відрізків. Також і двовимірна лінійна функція \(f(x,y) = px + qy + r\) інтегруватиметься точно за допомогою двовимірної формули прямокутників. Можна вибрати \(f(x,y) = 2 x + y\) і створити відповідну тестову функцію, яка автоматично перевірятиме наші дві реалізації двовимірної формули прямокутників. Для обчислення інтеграла від (f(x,y)) ми скористаємося SymPy, щоб уникнути помилок при ручному обчисленні. Тестова функція буде такою:
Як альтернативу методу прямокутників можна вибрати формулу трапецій. Висновок формули для подвійного інтеграла і реалізація використовує ті самі ідеї, які були наведені для формули прямокутників. Тому Вправа 3.13: Вивести та реалізувати формулу трапецій для подвійного інтеграла має бути виконанообов'язково.
Формула прямокутників для потрійного інтеграла
Так як нам вдалося узагальнити метод працюючий для одномірного інтеграла на подвійний, то, природно, досить просто поширити його на три виміри. Розглянемо потрійний інтеграл
Отримаємо формулу прямокутників цього інтеграла. Введемо допоміжні функції:
Для кожного з цих одновимірних інтегралів застосуємо формулу прямокутників:
Таким чином, отримаємо прямокутників для потрійного інтеграла:
Реалізація¶
Аналогічно випадку подвійного інтеграла, створимо файл rectangular_triple.py з реалізацією відповідних функцій:
Метод Монте-Карло длячисельного інтегрування з областей складної форми¶
Повторне використання формул для одномірного інтегрування при чисельному знаходженні подвійних і потрійних інтегралів можливе лише, якщо область інтегрування прямокутна (прямокутник або паралелепіпед). Для інших форм областей слід використовувати зовсім інший підхід. Найбільш загальний спосіб для двовимірних і тривимірних областей полягає в розбитті їх на велику кількість маленьких трикутників або тетраедрів та використання чисельного інтегрування для кожного трикутника або тетраедра. Повний алгоритм та його реалізація надто складні для даного спецкурсу. Натомість ми використовуватимемо альтернативний, досить простий і загальний метод, званий методом інтегрування Монте—Карло (ММК). Код його реалізації можна вмістити на півсторінки, але ця реалізація вимагає значно більше виконання функцій при обчисленні подвійних інтегралів у порівнянні з методом прямокутників.
Однак ММК є більш ефективним ніж формула прямокутників при обчисленні багатовимірних інтегралів (більше трьох змінних) за гіперкубічними областями. Подані ідеї для подвійних і потрійних інтегралів можуть бути легко узагальнені на чисельне інтегрування за змінними. І тут формула прямокутників вимагає обчислення \(m\) сум. При \(N\) часткових відрізках по кожній змінній формула прямокутників вимагатиме \(N^m\) обчислень функцій. Це означає, що обчислювальна робота зростає як експонента від кількості розмірностей. ММК не залежить таким чином від кількості розмірностей і є кращим при чисельному знаходженні багатовимірних інтегралів. Тому має сенс у розділі, присвяченій чисельному інтегруванню, розглянути застосування ММК як складних областей, так великихрозмірності.
Алгоритм інтегрування з ММК¶
Ідея ММК для чисельного інтегрування \(\int_a^b f(x) dx\) полягає у використанні теореми про середнє з математичного аналізу, в якій стверджується, що інтеграл \(\int_a^b f(x) dx\) дорівнює добутку довжини відрізка ( тут \(b - a\) ) і середнього значення \(\bar\) функції \(f\) на відрізку \([a, b]\). Середнє значення може бути обчислено за допомогою вибірки значень \(f\) на безлічі випадкових точок всередині області та обчисленні їх середнього арифметичного. У багатовимірному випадку, інтеграл оцінюється як добуток площі (обсягу) області та середнього значення функції, яке знову обчислюється за вибіркою на безлічі випадкових точок.
Введемо деякі величини, які дозволять нам формалізувати алгоритм чисельного інтегрування. Нехай дано двовимірний інтеграл:
де \(\Omega\) - двовимірна область задана за допомогою допоміжної функції \(g(x,y)\) :
Таким чином, межа області \(\partial\Omega\) задана неявною функцією (кривий) \(g(x,y) = 0\). Такий опис областей поширений останні десятиліття, у своїй \(g\) називається функцією рівня, а межа \(g = 0\) — нульовим контуром функції рівня. Для простих областей можна легко побудувати функцію (g) вручну, але в складніших промислових додатках слід звернутися до математичних моделей побудови (g).
Нехай \(A(Omega)) - площа області \(Omega) . Ми можемо чисельно знайти інтеграл у наступному ММК:
- поміщаємо область \(\Omega\) всередину прямокутника \(R\);
- генеруємо велику кількість випадкових точок на (R);
- обчислюємо частку \(q\) точок, які потрапили в область \(\Omega\);
- наближаємо \(A(Omega)/A(R)) числом q , тобто, вважаємо \(A(Omega) = qA(R)\);
- обчислюємо середнє значення \(\bar\) функції \(f\) на області \(\Omega\);
- обчислюємо наближене значення інтеграла як (A (Omega) bar).
Зазначимо, що площу \(A(R)\) прямокутника \(R\) легко обчислити, при тому що площа \(A(\Omega)\) нам не відома. Однак, якщо припустити, що частка площі \(A(R)\) займаною областю \(\Omega\) така ж як частка випадкових точок, що потрапили всередину \(\Omega\), можна отримати просте наближення для \(A(\) Omega)\).
Для того, щоб отримати уявлення про метод, розглянемо область, обмежену колом і поміщену прямокутник, як показано на малюнку нижче. Набір випадкових точок зображено синім кольором.
Реалізація¶
Функція на Python, що реалізує чисельне інтегрування \(\int_\Omega f(x, y) dxdy\) , може виглядати так (файл mc_double.py ):
Візьмемо найпростіший тестовий приклад: знайти площу прямокутника \([0, 2] \times [3, 4.5]\) , поміщеного у прямокутник \([0, 3] \times [2, 5]\) . Точне значення площі - 3, але ММК, на жаль, ніколи не є точним, тому неможливо передбачити результат алгоритму. Все, що ми знаємо - це те, що наближене значення інтеграла має прагнути до 3 при числі випадкових точок, що прагне до нескінченності. Також при фіксованому числі точок ми можемо запускати алгоритм кілька разів і отримувати різні значення, які коливаються близько до точного значення, так як різні вибірки точок використовуються при різних викликах алгоритму.
Площа прямокутника можна обчислити за допомогою інтегралу \(\int_0^2\int_3^ dxdy\) , тому в цьому випадку покладемо \(f(x,y) = 1\) , а функцію \(g\) можна визначити, наприклад, наступним чином: 1 якщо точка \((x, y)\) лежить усерединіпрямокутника \([0, 2] \times [3, 4.5]\) інакше -1. Нижче наведено приклад використання функції MonteCarlo_double в інтерактивній сесії IPython для обчислення площі з різними потужностями вибірки:
Ми бачимо, що значення коливаються близько 3, що підтверджує правильну реалізацію, але в принципі за неточною відповіддю можуть ховатися помилки.
Відомо, що стандартне відхилення наближеного значення інтеграла сходиться як \(n^\) , де \(n\) - потужність вибірки (кількість випадкових точок). Такий вид оцінки швидкості збіжності може використовуватися при тестуванні реалізації методу, але ми цей спосіб не розглядатимемо.
Тестова функція для функцій із випадковими числами¶
Щоб написати тестову функцію, нам потрібен юніт тест, який має ідентичну поведінку за різних запусків. Це складним під час використання випадкових чисел, оскільки ці числа змінюються щоразу під час запуску алгоритму, і кожен запуск дає різні результати. Стандартний спосіб тестування алгоритмів із випадковими числами полягає у фіксуванні початкового числа генератора випадкових чисел. Тоді послідовність чисел буде однією і тією самою при кожному запуску алгоритму. Припустимо, що функція MonteCarlo_double працює. Ми фіксуємо початкову точку, спостерігаємо деякий результат, і беремо цей результат як правильний. За умови, що тестова функція завжди використовує цю початкову точку, ми повинні отримувати точно такий же результат кожного разу, коли викликається функція MonteCarlo_double . Наша тестова функція може виглядати так (файл mc_double.py):
Інтеграл по колу¶
Тест наведений вище включав найпростішу функцію \(f(x,y) = 1\) . Нам слід виконати також тести для не постійної функції (f) і для більш складної області.Нехай \(\Omega\) - коло з центром на початку координат і радіусом 2 і нехай \(f(x,y) = \sqrt\) . Такий вибір дозволяє обчислити точне значення інтеграла: в полярних координатах \(\int_Omega f(x,y) dxdy\) перетворюється на \(2\pi\int_0^2 r^2dr = 16\pi/3\) . Ми повинні бути готові до досить грубих наближень, що коливаються біля цього результату. Ми знаємо з досвіду попереднього тестового прикладу, що точніші результати можна отримати на більшій кількості випадкових точок. Нижче наведено приклад тестової функції: