Дослідження точності чисельного диференціювання

 

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

У такий спосіб порядок точності результату відносно кроку сітки дорівнює числу залишених членів ряду, або, іншими словами, він дорівнює числу вузлів інтерполяції мінус порядок похідної. Тому мінімальне число вузлів необхідне для обчислення m -ої похідної, дорівнює m+1, воно забезпечує перший порядок точності.

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

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

 

Метод Рунге-Ромберга

 

Загальна ідея методу така: маємо деяку наближену формулу (х,к) для обчислення величини z (х) за її значеннями на рівномірній сітці з кроком h, а залишковий член цієї формули

.              (6.14)

Наприклад, , —задана функція. Нехай , ,

,

. Тут p = 2. Якщо скористатися тією самою наближеною формулою для обчислення значення z в точці х, але використовуючи сітку з кроком rh, дістанемо

(6.15)

Віднявши (6.14) від (6.15), дістанемо першу формулу Рунге для оцінки похибки

.        (6.16)

Перший доданок у (6.16) є головним членом похибки, тобто розрахунок на другій сітці дає змогу оцінити похибки на першій сітці з точністю до членів вищого порядку. Виключаючи за допомогою (6.16) величину   з (6.14), дістанемо другу формулу Рунге

,    (6.17)

яка дає результат з вищим порядком точності, ніж (6.14). Іноді уточненнярезультату за формулою (6.17) називають уточненням за Річардсоном. Розглянемо приклади застосування описаного вище процесу для підвищення точності в задачі чисельного диференціювання.

Приклад 1 Нехай функція   задана таблицею. Обчислити у' (3).

1 2 3 4 5	0,000 0,301 0,478 0,602 0,699

    Розв ’ язання. Скориставшись формулою  при , дістанемо . Збільшуючи крок удвічі (), дістанемо .

За формулою (6.16) при р = 2 , що лише на 2% відрізняється від шуканого значення у' (3)= 0,145.

Приклад 2 За допомогою методу Рунге вивести формулу чисельного диференціювання порядку  з формули більш низького порядку .

Розв ’ язання. Маємо

, .

Порядок точності цих формул , а коефіцієнт збільшення кроку , тому уточнення за методом Рунге дає формулу

.

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

Метод Рунге узагальнюється на довільну кількість сіток.

Приклад 3 За допомогою розвинення в ряд Тейлора для функції  і  дістаємо

. !

, .    (6.18)

Приклад 4 Для односторонньої різницевої похідної  при ,  маємо

, !, .

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

, , .

Це система лінійних рівнянь відносно величин  і , . Використавши формули Крамера, дістанемо уточнений розв'язок за формулою Ромберга

, (6.19)

де .

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

де .

Функції  мають, очевидно, такі дві властивості:

а) ,  - символ Кронекера;

б)  - дійсні коефіцієнти, тобто  є многочленами від . Тому функція

(6.20)

є інтерполюючою функцією для  (  - дійсні коефіцієнти), а величина

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

Оскільки система функцій   не при всіх  і не на довільному інтервалі буде системою Чебишева, то інтерполяційна функція (6.20)існує не для будь-якої послідовності . Але для послідовностей, які найчастіше трапляються на практиці, а саме:

а)  (послідовність Рунге - Ромберга);

б)  можна довести, що , і тим самим існування многочлена Р(х) гарантується.

Зауваження

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

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

,

і може бути обчислений в точці  за допомогою, наприклад, алгоритму Ньютона без розв'язування системи.

3 Якщо сітки такі, що , тобто згущення їх відбувається за одну і ту саму кількість разів, то зручніше застосувати рекурентно метод Рунге. Це робиться таким чином. Спочатку на кожній парі сіток  методом Рунге вилучають головний член похибки . Уточнені значення, таким чином, групуються в пари і далі вилучається похибка наступного порядку . Всього можна виконати  уточнень. При кожному уточненні за формулою (6.16) обчислюється апостеріорна оцінка точності.

4 Якщо формула для обчислення   має симетричний вигляд, то на рівномірній сітці часто всі непарні члени ряду (6.18) перетворюються на нуль (див. приклад 3). У такій ситуації користуватися формулою (6.19) невигідно. Потрібно залишити в сумі (6.18) члени   і відповідно змінити формулу Ромберга. Те саме стосується і рекурентної процедури Рунге - при черговому вилученні порядок точності підвищується на 2, а не на 1.

5 Число членів суми (6.18) пов'язане з кількістю неперервних похідних у функції, для якої обчислюються  і  (див. приклади 4, 5). Для не досить гладких функцій недоцільно брати велике число сіток для уточнення. Практично навіть для гладких функцій використовують не більше 3 – 5 сіток, причому, як правило, беруть відношення   кроків сіток, що дорівнює 2.

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

 

Процес Ейткена

Метод розрахунків на декількох сітках застосовується для підвищення порядку точності і в тому випадку, коли невідомий порядок головного члена похибки. Він має назву процесу Ейткена. Нехай

,

але  - невідоме. Проводяться обчислення на трьох сітках з кроками

.

Нехтуючи членами порядку , дістаємо

.

Звідси знаходимо . Далі можна скористатися вже відомим методом Рунге, який можна трактувати таким чином. Утворимо комбінацію  і виберемо  так, щоб

.

Дістанемо , причому

.

 

Питання і завдання до розділу 6

 

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

2 Апріорна оцінка похибки чисельного диференціювання.

3 Апостеріорна оцінка похибки чисельного диференціювання.

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

5     Скласти таблицю наближених значень похідної функції  за таблицею її значень

x	1	1.1	1.2	1.3
y	2	2. 144	2. 297	2. 462

6 Переконатися в тім, що формула чисельного диференціювання  має другий порядок точності.

7 Вивести формули чисельного диференціювання на основі лінійної інтерполяції.

8 Функція  задана таблично

х	1	2	3	4	5	6
у	1.34	3.05	5.21	7.22	9.11	12.34

Обчислити .

 

Р озділ 7

Чисельне інтегрування функцій

 

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

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

Рис. – 7.1

У практичних розрахунках нерідко виникає потреба в обчисленні визначених інтегралів вигляду

,

де функція та вагова функція  неперервні на відрізку .

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

Найчастіше підінтегральну функцію заміняють на деякий узагальнений поліном. Тоді внаслідок лінійності такої апроксимації функцію можна записати так:

,

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

,

де вузли; - ваги; похибка або залишковий член квадратурної формули.

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

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

Це відношення називається квадратурною формулою. При цьому: права частина виразу  називається квадратурною сумою. Тут параметри квадратурної формули: квадратурні (вагові) коефіцієнти; квадратурні вузли.

Якщо межі інтегрування  являються квадратурними вузлами, то отримуємо формулу замкненого типу. Інакше маємо квадратурну формулу відкритого типу.

Величина  називається похибкою квадратурної формули .

Якщо для деякої функції  маємо  то квадратурна формула являється для даної функції точною.

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

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

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

Лема. Степінь точності формули  не може бути вище за  при будь-якому виборі параметрів

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

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

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

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

Розглянемо формули для наближеного обчислення інтегралів

.                    (7.1)

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

Отже, нехай задана рівномірна сітка , , . Тобто крок  – величина постійна й розбиває відрізок  на  рівних інтервалів. Формули Ньютона-Котеса - формули замкненого типу. Позначимо . За наближену функцію  оберемо інтерполяційний поліном Лагранжа

 де .

Отже, заданий інтеграл може бути поданий у вигляді

Таку квадратурну формулу називають квадратурною формулою інтерполяційного типу.

Нехай  – виражена в сіткових кроках довжина . Тоді

;

.

У такому випадку ваги можна розрахувати так:

.                                        (7.2)

Формула (7.2) остаточно визначає ваги квадратурної формули Ньютона-Котеса. Замінимо в ній  і введемо позначення . Тоді коефіцієнти

                (7.3)

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

.                (7.4)

Для коефіцієнтів Котеса мають місце співвідношення:

1 .

2 .

З’ясуємо питання про степінь точності квадратурної формули .

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

Звідси можна зробити висновок, що квадратурні коефіцієнти  формули  є єдиним розв’язком лінійної системи рівнянь

яка отримана із  при

Розглянемо окремі випадки квадратурних формул Ньютона-Котеса з рівновіддаленими вузлами, в яких підінтегральна функція  замінена на інтерполяційний поліном Лагранжа різного степеня.

 

7.1.1 Формула середніх  (формула прямокутників)

 

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

().

Це найпростіша квадратурна формула. Розклавши  у ряд Тейлора довкола точки  

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

.

Формула середніх є точною для лінійної підінтегральної функції , оскільки тоді .

Природно, що точність формули для довільної  можна підвищити, якщо скористатися докладнішою сіткою

Це так звана складена формула середніх або формула прямокутників. У разі рівномірної сітки, тобто якщо , формула виглядатиме так:

Наведені оцінки R справедливі, якщо існує неперервна ; якщо ж  кусково-неперервна, має місце лише мажорантна оцінка

.

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

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

.

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

На докладнішій сітці  одержимо складену формулу трапецій:

На рівномірній сітці вона стає такою:

).

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

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

Квадратурна формула Симпсона є частковим випадком квадратурних формул Ньютона-Котеса при . Тут підінтегральна функція заміняється інтерполяційним поліномом Лагранжа 2-го степеня (рисунок 7.2). Із цієї причини формулу Симпсона ще називають формулою парабол.

Розіб’ємо відрізок  на 2 рівних відрізки і одержимо сітку , що містить три вузли. Формула Симпсона містить три коефіцієнти Котеса:

; .

.

; .

    (7.5)

Рис. – 7.2

Формула (7.5) є трьохточковою квадратурною формулою Симпсона.

Якщо  розбити на парну кількість відрізків, що дорівнює , і до кожного часткового здвоєного проміжку , ,…,  застосувати формулу Симпсона, то одержимо складену формулу Симпсона

.                     (7.6)

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

.

Отже, дійдемо формули парабол

.

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

Ця формула має такий залишковий член:

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

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

f(x):

//Повертає значення підінтегральної функції

End

//Повертає одну суму з формули Симпсона

//a – лівий кінець відрізка інтегрування

//h – крок

sum1(n,h,a):

1 temp:=0

2 for i:=1 to (an div 2) do

3 temp+=f(a+(2*k-1)*ah)

4 done

5 return temp

End

//Повертає одну суму з формули Симпсона

//a – лівий кінець відрізка інтегрування

//h – крок

sum2(n,h,a):

1 temp:=0;

2 for k:=2 to (an div 2) do

3 temp+=f(a+(2*k-2)*ah)

4 done

5 return temp

End

Calc_Integrate_Simpson(a,b,

1 n:=4

2 h:=(b-a)/n

3 repeat

4 I1:=(h*(f(a)+f(b)+4*sum1(n,h,a)+2*sum2(n,h,a)))/3;

5 n:=2*n;

6 h:=(b-a)/n;

7 I2:=(h*(f(a)+f(b)+4*sum1(n,h,a)+2*sum2(n,h,a)))/3

8      m:=abs(I1-I2)

9 until (m<eps)

10 return I2

End

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

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

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

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

       (7.7)

Отже, знайдемо параметри  з таких умов:

       (7.8)

Це система  нелінійних алгебраїчних рівнянь відносно , .

Для подальшого спрощення вважатимемо, що .

Якщо  одержимо  і система  набере вигляду

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

яка є точною для будь-якого полінома 1-го степеня.

Якщо , система  матиме такий вигляд ():

Розв’язавши цю систему, знайдемо:

тобто маємо квадратурну формулу

,

яка є точною для будь-якого полінома 3-го степеня.

За довільного  як вузли квадратурної формули Гауса беруть нулі поліномів Лежандра

а ваги цієї квадратурної формули визначають за таким виразом:

         (7.9)

Маючи значення  та вузлів  на відрізку , значення інтеграла на довільному відрізку  обчислюється за такою квадратурною формулою Гауса:

Похибка квадратурної формули Гауса має вигляд

Зауважимо, що, починаючи з , і вузли, і ваги є ірраціональними