Додаток 3

ОБЧИСЛЕННЯ ВИЗНАЧЕНИХ ІНТЕГРАЛІВ

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

Якщо функція f(x) безупинна на відрізку [а, b] і відома її первісна функція F(х), то визначений інтеграл цієї функції в межах від а до b може бути обчислений по формулі Ньютона-Лейбница

(1)

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

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

(2)

Сума називається квадратурною сумою, а рівності виду (2) являються формулами механічних квадратур. При їхній побудові можна розпорядитися 2n+1 параметром: n - абсцис x_i, n - коефіцієнтів A_i і останній параметр - число вузлів n. Усі ці параметри слід вибрати так, щоб формула (2) давала досить малу погрішність для широкого класу функцій f(x). Очевидно, що чим більше значення n, тим більшої точності можна досягти в рівності (2).

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

Виберемо на відрізку інтегрування [а, b] n + 1 довільних крапок х₀,х₁,..., х_n і інтерполюємо функцію f(x) по її значеннях у цих крапках за допомогою полінома Лагранжа

(3)

(4)

де L_i(X_i) = f(x_i), i =0,1,…,n... Точне значення інтеграла

(5)

Якщо інтерполяція (3) була досить точна, то залишковим членом

можна зневажити. Після цього одержуємо наближену рівність

(6)

де

(7)

Квадратурна формула (6) називається квадратурною формулою "закритого типу", якщо межі інтегрування а і b є вузлами інтерполяції, у противному випадку - "відкритого типу".

Помітимо, що:

коефіцієнти A_i при даному розташуванні вузлів не залежать від функції f(x);

для полінома ступеня n формула (6) - точна, тому що тоді Ln(x) f(x).

Квадратурні формули Ньютона-Котеса

Серед інтерполяційних формул раніше всіх були побудовані широко відомі і часто застосовувані квадратурні формули Ньютона-Котеса. Вони відносяться до випадку рівновіддалених вузлів інтерполяції.

Нехай для даної функції f(x) потрібно обчислити інтеграл . Вибравши крок h = (b - а) / n, розіб'ємо відрізок (а, b) за допомогою рівновіддалених крапок х₀= a, x_i= х₀ + ih, (i=1, ..., n - 1), x_n=b на n рівних частин.

Заміняючи функцію f(x) відповідним інтерполяційним поліномом Лагранжа, одержуємо наближену квадратурну формулу

(8)

де A_i - деякі постійні коефіцієнти.

Одержимо явні вирази для коефіцієнтів А_i формули (8). Позначимо

і

З огляду на те, що х_i = x₀ + ih, запишемо формулу (4) для інтерполяційного полінома Лагранжа в наступному виді:

(9)

Тому що q = (х - х₀) / h і dq = dx / h то, зробивши заміну перемінних у визначеному інтегралі (7), будемо мати:

(10)

Оскільки h = (b - a) / h, то звичайно думають

A_i = (b - a)H_i,

де

(11)

- постійними, що називаються коефіцієнтами Котеса.

Квадратурна формула (8) при цьому приймає вид

(12)

де х_i = а + ih, i = 0,1,..., n.

Неважко переконатися, що справедливо співвідношення

Формула трапецій

Застосовуючи формулу (11), при n = 1 маємо:

Звідси

(13)

Ми одержали відому формулу трапецій для наближеного обчислення визначеного інтеграла (рис. 1).

Рис. 1. Формула трапецій

Залишковий член (помилка) квадратурної формули (13) у припущенні, що f "(х) безупинна на відрізку [х₀, х₀ + h], пропорційний кроку в третьому ступені і виражається формулою

(14)

Оскільки відрізок інтегрування [а, b] може бути досить великим, то для обчислення інтеграла розділимо проміжок інтегрування на k рівних частин [x₀, x₁] ,[х₁, х₂],... ,[x_k-1,x_k] і до кожної з них застосуємо формулу трапецій. Думаючи h = (b - а) / (k _* n) і позначаючи через f_i=f(x_i) i=0,1,..., k значення функції f(x) в крапках х_i , будемо мати

(15)

Геометрично формула (15) виходить у результаті заміни графіка підінтегральної функції f(x) ламаною лінією.

Формула Симпсона

По формулі (11) при n = 2 одержуємо

Тому що х₂ - х₀ = 2h, маємо

(16)

Формула (16) має назву формули Симпсона. Геометрично ця формула виходить у результаті заміни; кривої f(x) параболою, що проходить через три крапки

A(х₀, f(x₀)) B(х₁, f(x₁)) і C(х₂, f(x₂)) (рис. 2).

Рис. 2. Формула Симпсона

Залишковий член формули Симпсона дорівнює

(17)

Для обчислення , як і у випадку формули трапецій, розіб'ємо проміжок інтегрування [а,b] на k рівних частин. При цьому відстань між вузлами інтегрування h = (b - а) / 2k. Застосовуючи формулу Симпсона до кожного під інтервалу довжиною 2h і позначаючи через f_i = f(x_i), i =0,1,…,k значення функції f(x) в крапках X_i, будемо мати

(18)

Звичайно довжина під інтервалу інтегрування 2h = (b - а) / k називається кроком інтегрування.

Квадратурна формула Ньютона

Виконавши обчислення для n = 3 по формулі (11), одержимо квадратурну формулу Ньютона

(19)

яка іноді називається "правилом трьох восьмих".

Для того щоб поширити правило "трьох восьмих" на інтеграл необхідно розділити відрізок інтегрування на k рівних частин. Тоді відстань між вузлами інтерполяції h = (b - а) / 3k. Застосовуючи до кожного побудованого відрізка [х₀, х₃], [х₃, х₆],..., [х_3k-3, х_3k] формулу "трьох восьмих", одержимо

,

де f_i = f(x_i),х_о=а, x_i=х₀+ih, i=0,1,...,3k.

Формули Ньютона-Котеса з великим числом вузлів порівняно рідко застосовуються в обчислювальній практиці. Звичайно воліють користатися формулами з малим значенням n. Для зменшення погрішності результату попередньо розбивають відрізок [а, b] на досить малі інтервали і до кожного з них застосовують квадратурну формулу з малим числом вузлів.

Приведемо обчислені по формулі (11) коефіцієнти Котеса для n від 1 до 8:

n = 1: Н₀ = H₁ = 0.5;

n = 2: Н₀ = Н₂ = 0.1666667, H₁ = 0.6666667;

n = 3: Н₀ = Н₃ = 0.125, Н₁ = Н₂ = 0.375;

n = 4: Н₀ = Н₄ = 0.0777778, Н₁ = Н₃ = 0.3555556, Н₂ = 0.1333333;

n = 5: Н₀ = Н₅ = 0.0659722, H₁=Н₄ = 0.2604166, Н₂ = Н₃ = 0.1736111;

n = 6: H₀ = H₆= 0.0488095, H₁ = H₅ = 0.2571428, H₂ = H₄ = 0.0321428,

H₃ = 0.3238095;

n = 7: H₀ = H₇ = 0.0434606, H₁ = H₆ = 0.2070023, H₂ = H₅ = 0.0765626,

H₃ = H₄ = 0.1729745;

n = 8: H₀ = H₈ = 0.0348853, H₁= H₇ = 0.2076895, H₂ = H₆ = -0.0327336,

H₃ = H₅ = 0.3702292, H₄ = -0.1601410.

Квадратурна формула Гауса

Загальна квадратурна формула

(21)

при фіксованому числі вузлів n містить 2n параметрів A_i і x_i (i = 1,..., n). Розглянуті квадратурні формули

Ньютона-Котеса були отримані для довільно розташованих вузлів і точні для цілих функцій, ступінь яких не вище n-1. Для збільшення ступеня точності можна ще скористатися n вузлами x_i. Оптимальний вибір вузлів інтерполяції дозволяє збільшити ступінь точності формули (21) на n одиниць і зробить її вірної для багаточленів ступеня ≤2n – 1. Ця задача на основі використання інтерполяційних поліномів Лежандра була уперше вирішена Гаусом. Опускаючи висновок, приведемо загальну формулу Гауса для обчислення визначеного інтеграла

(22)

де x_i=a+(b-a)x_i, i=1,2,…,n...

Нижче приведені абсциси х_i і коефіцієнти A_i формули Гаусса для значень меж інтегрування а = 0, b = 1, при яких формула Гауса має вид:

(23)

Абсциси і коефіцієнти формули Гаусса

n = 1: X₁ = 0.5; А₁ = 1.0;

n = 2: X₁ = 0.2113249, x₂ = 0.7886751; А₁ = А₂ = 0.5;

n = 3: X₁ = 0.1127017, x₂ = 0.5, x₃ = 0.8872983;

A₁ = A₃ = 0.2777778, A₂ = 0.4444444;

n = 4: X₁ = 0.0694318, X₂ = 0.3300095, X₃ = 0.6699905,

A₁ = A₄ = 0.1739274, A₂ = A₃ = 0.3260726;

n = 5: X₁ = 0.0469101, X₂ = 0.2307653, X₃ = 0.5, X₄ = 0.7692347,

X₅ = 0.9530899;

А₁ = A₅ = 0.1184634, A₂ = A₄ = 0.2393143, A₃ = 0.2844444;

n = 6: X₁ = 0.0337652, X₂ = 0.1693953, X₃ = 0.3806904,

X₄ = 0.6193096, X₅ = 0.8306047, X₆ = 0.9662348;

A₁ = A₆ = 0.0856622, A₂ = A₅ = 0.1803808,

A₃ = A₄ = 0.2339570;

n = 7: X₁ = 0.0254460, X₂ = 0.1292344, X₃ = 0.2970774, X₄ = 0.5;

X₅ = 0.7029226, X₆ = 0.8707656, X₇, = 0.9745540;

A₁ = A₇ = 0.0647425, A₂ = A₆ = 0.1398527,

A₃ = A₅ = 0.1909150, A₄ = 0.2089796;

n = 8 X₁ = 0.0198551, X₂ - 0.1016668, X₃ = 0.2372338,

X₄ = 0.4062827, X₅ = 0.5917173, X₄ = 0.7627662,

X₇ = 0.8983332, X₈ = 0.9801449;

A₁ = A₈ = 0.0506143, A₂ = A₇ - 0,1111905,

A₃ = A₆ = 0.156853Е, A₄ = A₅ = 0.1813419.

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

Так, розбивши відрізок інтегрування [а, b] на k рівних частин і застосувавши до кожної частини формулу Гауса для фіксованого n, одержимо

(24)

де h = (b - а) / k - крок інтегрування;

- i-а абсциса j-го кроку інтегрування:

х₀ = а, х_j = х_j-1+ h, j = 1,2,..., k;

= x_j-1+x_ih, i = 1,2,…, n;

X_i, A_i - абсциси і коефіцієнти формули Гауса.

При однаковому числі ординат формула Гауса дає велику точність, ніж інші формули, але вона менш ефективна в методах з автоматичним вибором кроку інтегрування. Це порозумівається тим, що процедура автоматичного вибору кроку інтегрування заснована на дворазовому обчисленні інтеграла на кожнім кроці: спочатку з кроком h, а потім із кроком h/2. Отримані результати порівнюються і, якщо Δs= |s₁- s₂| менше заданої наперед припустимої помилки, то вважається, що інтегрування на даному кроці виконано правильно. У противному випадку крок поділяється навпіл і повторюється описана вище процедура інтегрування. Оскільки в квадратурних формулах Ньютона- Котеса відстань між вузлами інтегрування однакова, ординати, обчислені при кроці h, можуть бути використані і при половинному кроці. У формулах Гауса цього зробити не можна, тому що відстані між абсцисами не однакові. Методи автоматичного вибору кроку інтегрування (адаптивні квадратурні програми) будуть розглянуті нижче.

Найпростішою квадратурною формулою є формула Гауса при n = 1:

Х₁ = 0.5, A₁ = 1.0, що називається формулою прямокутників. Відповідно до цієї формули

(25)

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

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

(26)

де h = (b - а) / 2 ; х₁ = а + h / 2, х_j = х_j-1 + h, j = 2,3,...,k.

Рис.3. Формула прямокутників

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

Рис.4. Формула "лівих" і "правих" прямокутників

Ці формули звичайно називаються формулами "лівих" і "правих" прямокутників і відрізняються від прямокутників (26) тим, що у формулі "лівих" прямокутників

Х₁ = а, а у формулі "правих" прямокутників х₁ = а + h.

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

Адаптивні квадратурні програми

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

Користувач подібної програми вказує інтервал [а, b], заготовляє функцію, що обчислює f(x) для будь-якого х з цього інтервалу, і вибирає припустиму точність ε. Програма намагається обчислити величину Q таку, що

У процесі обчислень інтервал [а, b] розбивається на під інтервали [x_i,x_i+1]. У більшості програм кожен під інтервал виходить розподілом навпіл під інтервалу, отриманого на більш ранньому етапі обчислень. Реальне число під інтервалів, так само як їхнє розташування і довжини, залежать від функції f(x) і необхідної точності ε. Нижче довжина під інтервалу буде позначатися через h_i = x_i+1 – х_i.

Типовим прийомом є застосування до кожного під інтервалу двох різних квадратурних формул. Позначимо відповідні два результати через p₁ і Q₁ Наприклад, схеми, засновані на формулі Симпсона, використовують основну формулу (два елементарних відрізки)

і складену (чотири елементарних відрізки)

Вираження р_i і Q_i є наближеннями до інтеграла

Основна ідея адаптивної квадратури складається в порівнянні двох наближень P₁ і Q_i і одержанні при цьому оцінки їхньої точності. Якщо точність прийнятна, то один з результатів (краще Q_i) приймається за значення інтеграла на даному під інтервалі. Якщо точність недостатня, то під інтервал поділяється навпіл.

Типовою оцінкою точності є перевірка нерівності

де P - показник ступеня похідної, що входить у залишковий член квадратурної формули. Так, для формули Симпсона згідно (17) р – 1 = 4; h_i - поточний крок інтегрування; ε – задана користувачем припустима абсолютна точність обчислень, тобто

Приклади рішення задач

Задача А

Скласти функції обчислення визначеного інтеграла по формулі Ньютона-Котеса і Гауса для n=4. Одержати і порівняти методичні помилки наближеного обчислення інтеграла

по формулах Ньютона-Котеса і Гауса. Досліджувати вплив кроку інтегрування на точність обчислень.

Первісна функція заданого інтеграла дорівнює

Лістінг програми

#include <stdio.h>

#include <math.h>

#include <conio.h>

typedef double D;

int kf;

//функція fun() – підінтегральна функція

D fun(D x)

{

D cx;

kf ++ ;

сх=соs(х);

return cx*cx/sin(x) ;

}

// функція pf () - первісна функція

D px (D x)

{return cos(x)+log(tan(x/2));}

// функція n_k_4() – формула Ньютона -Котеса для n=4

D n_k_4(D a, D b, D f (D), int k)

{

D h, hk, s, xi;

Static D H0=0. 07777778,

H[4]={ 0.3555556, 0.1333333, 0.3555556, 0.1555555 };

int i, im=4*k;

hk=(b-a)/k;

h=hk/4 ;

s=(f (a)+f (b))*HO;

for ( i=1; i<im; i++)

{

xi=a+1*h;

s+=H[(i-1)%4]*f(xi);

}

return hk*s;

}

// функція gauss_4() - формула Гауса для n=4

D gauss_4 (D a, D b, D f(D), int k)

{

D h, xi, x, s;

int i, j;

static D A[4]={ 0.1739274, 0.3260726, 0.3260726, 0.1739274 },

X [4]={ 0.0694318, 0.3300095, 0.6699905, 0.9305682 };

h=(b-a)/k;

for(s=i=0; i<k; i++)

{

xi=a+i*h;

for(j=0; j<4; j++)

{

x=xi+X[ j]*h;

s+=f(x)*A[j];

}

}

return h*s;

{

void main()

{

D a=0.2, b=2.1, st, sp_n, ds_n, sp_g, ds_g;

Int k, kf_n,kf_g;

clrscr( );

printf("\n\t\t Залежність помилки обчислень від числа кроків" );

printf("\n k\t st\t kf_n\t sp_n\t ds_n\tkf_g sp_g\t ds_g");

st=pf(b) - pf(a);

for(k=l; k<9; k++)

{

kf =0;

sp_n=n_k_4(a, b, fun, k);

kf_n=kf;

ds_n=sp_n-st;

kf =0 ;

sp_g=gauss_4 (a, b, fun, k);

kf_g-kf;

ds g-sp g-st;

printf("\n %2d %11.3e %4d %11.3e %11.3e %5d %11.3e",

k, st, kf_n, sp_n, ds_n, kf_g, sp_g, ds_g);

}

getch();

}

Залежність помилки обчислень від числа кроків:

K	St	kf_n	sp_n	ds_n	kf_g	sp_g	ds_g
1	1.370е+00	5	1.465е+00	9.520е-02	4	1.352e+00	-1.816з-02
2	1.370е+00	9	1.385е+00	1.441е-02	8	1.368е+00	-1.958e-03
3	1.370e+00	13	1.374е+00	3.773e-03	12	1.370е+00	-3.783e-04
4	1.370е+00	17	1.371е+00	1.294e-03	16	1,370е+00	-9.969е-05
5	1.370е+00	21	1.371е+00	5.262e-04	20	1.370е+00	-3.211е-05
6	1.370e+00	25	1.370е+00	2.409е-04	24	1.370е+00	-1.194e -05
7	1.370е+00	29	1.370е+00	1.205е-04	28	1.370e+00	-4.946е-06
8	1.370е+00	33	1. 370е+00	6.469е-05	32	1.370е+00	-2.233e-06

Вище позначено:

k - число кроків інтегрування;

st - точне значення інтеграла;

kf_n, kf_g - число звертань до підінтегральній функції f(x) у методах Ньютона-Котеса і Гауса відповідно;

sp_n, sp_g - наближене значення інтеграла, отримане по формулах Ньютона- Котсса і Гауса;

ds_n, ds_g - помилка обчислення інтеграла в методах Ньютона-Котеса і Гауса.

Завдання на програмування

Задача 1 Скласти функції обчислення визначеного інтеграла по заданих квадратурних формулах Ньютона-Котеса і Гауса.

Використовувати складені функції при обчисленні заданих інтегралів. Досліджувати вплив кроку інтегрування на точність обчислень інтеграла

Для методів Ньютона-Котеса і Гауса діапазон зміни числа кроків однаковий і дорівнює k = кн(Δk)k_k, де кн - початкове, k_к - кінцеве значення числа кроків, Δk - крок зміни числа кроків інтегрування.

Діапазон зміни числа кроків інтегрування є функцією кількості вузлів:

n кн ( k) k_k

1 18 (2) 28

2 10 (2) 20

3 8 (2) 16

4 6 (2) 16

5,6 6 (1) 12

7,8 4 (1) 8

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

Результати обчислень повинні містити:

- число кроків інтегрування;

- число звертань до функції обчислення підінтегральній функції для методів Ньютона-Котеса і Гауса;

- точне значення інтеграла;

- наближені значення і помилки обчислення інтеграла по формулах Ньютона-Котеса і Гауса.

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

Варіанти визначених інтегралів, що обчислюються, і методи обчислень представлені в табл. 1 і 2 відповідно (Додаток 4).