Квадратурна формула Чебишева

Візьмемо за основу формулу і будемо вважати всі квадратурні коефіцієнти однаковими: . Тоді

. (7.10)

Параметри виберемо так, щоб формула була точною для всіх поліномів степеня не вище за При цьому достатньо розглянути функції

Для маємо

Покладаючи в приходимо до системи нелінійних рівнянь для визначення квадратурних вузлів

Отримана квадратурна формула називається формулою Чебишева. Зокрема при

7.3 Стійкість квадратурного процесу. Оцінки похибки

Стійкість квадратурних формул характеризує їх чутливість до різного роду похибок. Вона безпосередньо пов'язана з поняттям збіжності квадратурних формул.

Квадратурна формула буде збіжною за умови, що залишок при .

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

. (7.11)

Це неусувна похибка квадратурної формули.

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

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

Існує можливість оцінити похибку квадратурної формули до початку розв’язання задачі. Така оцінка називається апріорною. Оцінка похибки після розв’язання задачі називається апостеріорною.

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

.

Оскільки функція , то:

;

;

,;

,;

,;

.

Залишковий член загальної формули Симпсона

. (7.12)

Оскільки – неперервна на функція, то знайдеться така точка , що

.

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

У зв'язку з цим широкого застосування набуло правило Рунге апостеріорної оцінки похибки, суть якого полягає в тому, щоб, організувавши обчислення двох значень інтеграла на двох множинах вузлів, їх порівняти й одержати оцінку похибки. Найпоширеніше обчислення інтеграла двічі - із кроками та .

Якщо – точне значення інтеграла, – його наближене значення, обчислене з кроком , а – наближене значення інтеграла, обчислене із кроком , то похибки кожної квадратурної формули із кроком і можна записати відповідно у вигляді , , де – порядок точності формул. Обчислимо наближене значення інтеграла за однією квадратурною формулою спочатку із кроком , а потім із кроком . Одержимо

, .

; .

Одержали оцінку похибки методом Рунге

. (7.13)

Користуючись цією формулою можна уточнити наближене значення інтеграла

. (7.14)

Цю формулу називають формулою екстраполяції за Річардсоном.

Для формули Симпсона .

7.4 Вибір квадратурних формул

чисельного інтегрування

Ми одержали ряд формул чисельного інтегрування. Виникають запитання: яку формулу потрібно застосовувати в тому або іншому випадку, які формули більш вигідні і які менш вигідні. На ці питання не можна відповісти однозначно. Усе залежить від того, яким способом задана підінтегральна функція, які обчислювальні засоби використовуються, яка необхідна точність і таке інше.

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

Якщо обчислення ведуться вручну або за допомогою малих обчислювальних машин, то мають значення формули, що містять різниці. Менш вживані формули Гауса і Чебишева, тому що обчислення з багатозначними коефіцієнтами й абсцисами в цьому випадку є громіздкими. З формул, що не містять різниці, найчастіше застосовується формула Симпсона.

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

Тут необхідно зробити деякі зауваження щодо більш точних і менш точних формул. Ці терміни були введені нами при виведенні формул чисельного інтегрування й у них вкладався певний зміст. Потрібно чітко розуміти, що більш точна в цьому змісті формула не завжди дає практично більш точний результат. Справді, візьмемо найбільш точну з формул - формулу Гауса . Вона має вигляд

, (7.15)

де коефіцієнти С_i і абсциси x_i фіксовані і залежать тільки від n і [а,b]. Може трапитися, що підінтегральна функція набуває нульового значення у кожній із точок х_i, а абсолютна величина інтеграла від неї велика. Тоді різниця між точним значенням інтеграла і наближеним, отриманим за формулою Гауса, буде також дуже велика. У зв'язку з цим потрібно відзначити, що при виборі тієї або іншої формули чисельного інтегрування буває доцільно вивчити поводження підінтегральної функції і порівняти його з поводженням інтерполяційного многочлена, інтегруванням якого і знаходиться формула чисельного інтегрування. Іноді виникає необхідність розбивати відрізок інтегрування на окремі ділянки так, щоб краще описати поводження функції інтерполяційними многочленами.

Стосовно вибору квадратурних формул, то, очевидно, що за існування четвертої неперервної похідної від підінтегральної функції краще користуватися формулою парабол, за існування лише другої похідної та аналітичного задання - формулою середніх, а за табличного–формулою трапеції. Для функції високої гладкості найзручнішою є формула Гауса.

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

Швидкість збіжності квадратурної формули визначається оцінкою залишкового члена. Тобто, якщо , то квадратурна формула називається збіжною з р – порядком збіжності (нагадаємо, що для формули середніх, формули трапеції р=2, а для формули парабол р=4). Як правило, підінтегральна функція при цьому повинна мати похідну, яка входить до залишкового члена. Доведено, що чисельне інтегрування стійке за вхідними даними, хоча квадратурні формули нестійкі відносно похибок округлення, але ця нестійкість слабка і виявляється лише за розрахунків з малою кількістю правильних цифр.

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

х			Трапецій	Симпсона	Ейткена
0,00 0,25 0,50 0,75 1,00	0,0000 0,5000 0,7071 0,8660 1,0000	1,00 0,50 0,25	0,5000 0,6036 0,6433	–0,6381 0,6565	– – 0,6446

Як бачимо, дві формули дають результати невисокої точності. Погана точність формули Симпсона означає, що формула трапецій має фактично не другий порядок точності й уточнення методом Рунге тут недоцільне. А уточнення першого стовпця таблиці процесом Ейткена значно покращує результат; водночас з’ясовується, що в даному прикладі ефективний порядок точності формули трапецій .

Ефективний порядок точності виявився не цілим числом! З цим доводиться стикатися, якщо функція має особливості, а формула інтегрування явно цього не враховує, або якщо особливості має сама формула (це можливо в нелінійних формулах інтегрування).

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

Для обчислень можна використати пакет Maple:

> F1:=0.5;

> F2:=0.6036;

> F3:=0.6433;

> F:=F1+(F1-F2)^2/(2*F2-F1-F3);

> q:=0.5;

> beta=F-F1;

> p=ln(q)^(-1)*ln((F3-F2)/(F2-F1));

.

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

Використавши один із методів оцінки похибки квадратурних формул, можна оцінити точність проведення обчислень для поставленого завдання, якщо необхідно – вибрати оптимальний варіант обчислень, тобто найбільш точну формулу для даного випадку.

7.5 Чисельне інтегрування кратних інтегралів

 

Розглянемо K-вимірний інтеграл вигляду

(7.16)

де - деяка K-вимірна точка. Далі розглянемо подвійні інтеграли (K=2), оскільки їх можна інтерпретувати графічно.

Кубатурні формули, або формули чисельних кубатур, призначені для чисельного визначення кратних інтегралів.

 

 

Нехай функція визначена й неперервна в деякій обмеженій області . У цій області вибирається система точок (вузлів) . Для обчислення інтеграла приблизно покладемо

(7.17)

Щоб знайти коефіцієнти , зажадаємо точного виконання кубатурної формули (7.17) для всіх поліномів

(7.18)

степінь яких не перевищує заданого числа . Для цього необхідно й достатньо, щоб формула (7.17) була точною для добутку степенів

. Покладаючи в (7.17) , маємо:

(7.19)

Таким чином, коефіцієнти формули (7.17) можуть бути визначені із системи лінійних рівнянь (7.19).

Для того щоб система (7.19) мала розв’язок, необхідно, щоб число невідомих дорівнювало числу рівнянь. У випадку подвійного інтеграла () одержуємо

7.6 Вибір кубатурних формул

 

Для одержання заданої точності у разі К-кратного інтеграла сітковим (різницевим) методом потрібно виконати близько обчислень підінтегральної функції, де р – порядок точності сіткової формули.

Отже, якщо , вигідні сіткові методи, якщо ж то вигідний метод Монте-Карло. Так, за р=2 тривимірний інтеграл обчислюють сітковими методами, а при K=5 – методом Монте-Карло.

Розглянемо інтеграл по k-вимірній області, яка розбита сіткою на комірки (Рис. 7.4). Його можна обчислити послідовним інтегруванням:

Кожний однократний інтеграл легко обчислюється на даній сітці за квадратурними формулами типу

Послідовне інтегрування в усіх напрямках приводить до кубатурних формул, які є прямим добутком одновимірних квадратурних формул:

(7.20)

Наприклад, при k=2, якщо по кожному напрямку обрана узагальнена формула трапецій, а сітка рівномірна, то ваги кубатурної формули дорівнюють відповідно для внутрішніх, граничних і кутових вузлів сітки. Легко показати, що для двічі неперервно диференційованих функцій ця формула має другий порядок точності, і до неї застосуємо метод Рунге-Ромберга.

Взагалі для різних напрямків можна використати квадратурні формули різних порядків точності . Тоді головний член похибки має вигляд

Бажано для всіх напрямків використовувати квадратурні формули однакового порядку точності.

Можна підібрати ваги й положення ліній сітки так, щоб одновимірна квадратурна формула була точною для многочлена максимального степеня, тобто була б формулою Гауса. Тоді для випадку k=2:

(7.21)

де -нулі многочленів Лежандра й відповідні ваги. Ці формули розраховані на функції високої гладкості й дають для них більшу економію за кількістю вузлів у порівнянні з простішими формулами.

 

Метод послідовного інтегрування можна застосовувати до області довільної форми, наприклад, із криволінійною границею. Розглянемо цей випадок при K=2. Для цього проведемо через область хорди, паралельні осі , і на них уведемо вузли, розміщені на кожній хорді так, як нам потрібно (рис. 7.5). Представимо інтеграл у вигляді

Спочатку обчислимо інтеграл по уздовж кожної хорди за будь-якою одномірною квадратурною формулою, використовуючи введені вузли. Потім обчислимо інтеграл по ; тут вузлами будуть служити проекції хорд на вісь ординат.

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

Тоді друге інтегрування виконується за формулами Гауса

(7.22)

де , а й -нулі й ваги многочленів Чебишева другого роду.