Точність і стійкість чисельних методів
При чисельному розв’язуванні задач методом Ейлера єдиним параметром, яким можна управляти, є крок h. Типові програми потребують щоб користувач задав допуск EMBED Equation.3 . Ця величина використовується при виборі кроку інтегрування. Вибравши крок h достатньо малим, ми надіємось зменшити похибку до величини, яка не перевищує EMBED Equation.3 . Щоб це здійснити нам необхідно вміти оцінювати точність обчислень.
Запишемо формулу Тейлора, поклавши t=tk і h=hk :
EMBED Equation.3
Як правило, EMBED Equation.3 . Єдиним виключенням є вузол t0 , де y0=y(t0). Якщо з цього виразу виключити формулу методу Ейлера, отримаємо
EMBED Equation.3
У лівій частині записана похибка у вузлі tk+1 . Вона називається глобальною похибкою; саме цю величину нам необхідно контролювати. Глобальна похибка дорівнює різниці між обчисленим і “теоретичним” розв’язком. Є два джерела цієї похибки.
(1)Локальна похибка. Якщо по випадковості EMBED Equation.3 , то в правій частині нашого виразу дві різниці перетворюються в нуль. Але глобальна похибка нулю не дорівнює. Член, який залишився називається локальною похибкою або похибкою відсікання. Він в точності співпадає з залишковим членом EMBED Equation.3 у формулі Тейлора. Локальна похибка – це похибка, яка виникає при переході від tk до tk+1 , вважаючи, що всі величини при tk (і при t< tk) обчислені точно; в цьому понятті вона є штучною характеристикою.
(2)Поширювана похибка. Згідно теореми про середнє, другий доданок в правій частині можна записати у вигляді
EMBED Equation.3 .
Запис EMBED Equation.3 означає, що якобіан береться в деякій, наперед невідомій, точці. Отже, (1) і (2) можна сумувати наступним чином:
глобальна похибка при tk+1=(1+hkJ) (глобальна похибка при tk)+ локальна похибка при tk+1 .
Глобальна похибка у вузлі tk, y(tk)-yEMk множиться на величину (1-hkJ), яка називається множником переходу. Якщо EMBED Equation.3 , то похибка в методі Ейлера не збільшується, а якщо цей модуль більший одиниці, то похибки будуть збільшуватися. В останньому випадку говорять, що метод Ейлера нестійкий.
Нерівність EMBED Equation.3 еквівалентна нерівностям –2<hJ<0, які задають інтервал стійкості метода Ейлера. Якщо однорідне диференціальне рівняння (ОДР) нестійке, то ця нерівність ніколи не виконується і нестійкість методу Ейлера завжди виявляється. Якщо ж ОДР стійке (J<0), то ці дві нерівності ведуть до умови h<-2/J. Але навіть при стійкості задачі метод Ейлера може виявитися нестійким, якщо крок інтегрування не вибраний меншим ніж EMBED Equation.3 .
Приклад. Стійкість метода Ейлера.
EMBED Equation.3
При t<1 рівняння нестійке, що веде до збільшення похибки методу Ейлера, оскільки EMBED Equation.3 при будь-якому h. При t>1 рівняння є стійким. Але, не дивлячись на стійкість рівняння, метод Ейлера буде нестійким, якщо не виконано обмеження на крок EMBED Equation.3 . Прийнявши EMBED Equation.3 і h=0.1, ми повинні чекати збільшення похибки при t>3. На рис.1 точки (tk, yk) з’єднанні ламаними. Нестійкість спочатку “набирає розмаху” і вже з точки t=4 проявляється дуже різко. Для даної задачі множник переходу від’ємний, тому глобальна похибка міняє знак від вузла до вузла. Подібний ефект майже завжди є наслідком нестійкості чисельного метода. В таких випадках не потрібно зв’язувати отриманий чисельний розв’язок з досліджуваним фізичним процесом.
Дуже важливо розуміти, що чисельний метод може бути нестійким, навіть якщо ОДР стійке. Є задачі, розв’язок яких мало змінюється для всіх t>0. Якщо для задачі EMBED Equation.3 , EMBED Equation.3 =1000, то для методу Ейлера обмеження стійкості на крок інтегрування має вигляд h<2 10-3. Якщо нас цікавить розв’язок в пограничному шарі, то такий малий крок має фізичний зміст але в області плавної зміни розв’язку вибір такого кроку не має змісту. В першому випадку говорять, що крок інтегрування для чисельного методу визначається із умови точності, а в 0
100
200
100
0
1
2
3
4
5
Рис.1. Метод Ейлера: нестійкість.
Рис.2. Метод Ейлера: поведінка при розв’язку жорсткої задачі.
6
-2
-10
14
0
0.4
0.8
1.2
1.6


10
2
-6
2.4
2.8
другому випадку – із умови стійкості. Можна побачити, які проблеми виникають при використанні методу Ейлера для розв’язку жорстких задач, якщо вивчити рис.2. Чим дальше ми відхиляємося від розв’язку, тим більше “обчислений” нахил відрізняється від правильного значення. Для розв’язку жорстких задач метод Ейлера використовується рідко, якщо нас не цікавить в першу чергу пограничний шар.
Крок h входить як в локальну похибку h2y”, так і в множник переходу 1+hJ. Вибравши h достатньо малим, ми можемо зробити локальну похибку скільки потрібно малою і у випадку стійкості задачі досягнути того, щоб множник переходу був меншим одиниці. Але ні y”, ні J в явному вигляді не відомі. Більшість програм стараються оцінити y”, але ігнорують J. Інколи говорять, що такі програми здійснюють тільки контроль локальної похибки.
Процедура контролю похибки має на меті вибрати крок приблизно обернено пропорційний квадратному кореню з величини другої похідної, тому вона здійснює операції з локальною похибкою. Але локальна похибка – це лише частина загальної картини; на величину глобальної похибки суттєво впливає множник переходу. Знаючи точний розв’язок задачі, ми можемо обчислити у”(t) безпосередньо. Для EMBED Equation.3 функція зменшується по t, а при EMBED Equation.3 функція ніколи не перевищує 1. Таким чином, стратегія вибору кроку, яка здійснює операції з цими функціями, веде до послідовного зростаючого h для першого випадку і більш менш постійного h для другого випадку. В двох цих випадках виникає протиріччя з умовою стійкості для методу Ейлера, результатом якого буде швидка і чітка нестійкість. Але якщо ми звернемося до алгоритму вибору кроку, який не використовує другу похідну явно, а використовує тільки її оцінку, то переконаємося, що такий алгоритм не породжує “вибухів” і гігантських похибок. Навпаки, h підтримується приблизно пропорційним EMBED Equation.3 для двох випадків, а похибка контролюється. Підкреслимо, що крок буде малим навіть в області плавної зміни розв’язку жорсткої задачі. Це залежить від властивостей сімейства розв’язків.
Оцінка другої похідної базується на обчисленні різниці перших похідних y’n-1 і y’n. Якщо ми знаходимося на кривій розв’язку, то ці криві точні і оцінка другої похідної мала, бо всі y’(t) малі. Але EMBED Equation.3 і EMBED Equation.3 . Тому ці два числа необхідно розглядати як точні значення похідних, але від різних розв’язків. Але у випадку жорсткої задачі найбільше збурення кривої розв’язку може породити сильну зміну похідної y’ (оскільки J велике). Отже, різниця між значеннями повинна бути велика, а це веде до малого кроку інтегрування.
Приклад. Змінний крок метода Ейлера для жорсткої задачі.
Розглянемо жорстку задачу
EMBED Equation.3
Виходячи з точного початкового значення при t=0.5, попробуємо проінтегрувати рівняння до t=1.0 методом Ейлера з простою стратегією вибору кроку. Візьмемо початковий крок рівним 0.05 і задамо допуск на похибку EMBED Equation.3 . Якщо крок від t до t+h виявився невдалим, тобто оцінка похибки була велика, то попробуємо ще раз з кроком h/2. Цей алгоритм легко запрограмувати. Для досягнення точки t=1 достатньо 669 обчислень функції f(t,y) , і відносна похибка в останній точці рівна 0.005. Отже, середній крок рівний 1/(2669), або приблизно 0.0007; іншими словами, методу Ейлера з постійним кроком h=0.0007 потрібна була б та ж кількість обчислень функції f. Даний крок значно менший від того, який, за інтуїтивним відчуттям, потрібний, щоб дослідити повільно змінну функцію, і це пояснюється жорсткістю рівняння.
Поняття стійкості чисельного метода відіграє фундаментальну роль в наших дослідженнях. Програми з вибором постійного кроку можуть породжувати нестійкість, що видно з прикладу 1. При виборі змінного кроку нестійкість виникає порівняно рідко, але такі програми можуть бути неефективними, на що вказує вище приведений приклад.
Резюме. ОДР є стійкими чи нестійкими незалежно від того, яким чином воно розв’язується. Чисельний метод може бути стійким чи нестійким для конкретної задачі в залежності від розміру кроку.