Головна
1. Метод дихотомії
2. Метод простої ітерації
3. Метод Ньютона
4. Метод Ньютона для кратних коренів
5. Обчислення коренів поліномів методом Ньютона
6. Визначення екстремальних точок функції методом Ньютона
7. Точність обчислень у разі застосування методу Ньютона
8. Метод Січних
9. Метод хорд
10. Метод Мюллера
Розділ 4. Роз'язування нелінійних рівнянь

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

f(x) = 0     (4.1)

виконується переважно чисельними методами, основними властивостями яких є ітераційність рішення і локальність апроксимації.
   На відміну від лінійних рівнянь не існує прямих методів розв'язання нелінійних рівнянь. Загалом процедура розв'язання нелінійних рівнянь складається з двох етапів: попереднє знаходження інтервалів, що містять лише один корінь (локалізація коренів) і подальше уточнення коренів (розв'язання рівняння).
   Процедура розв'язання починається з вибору початкової точки хо і обчислення нев'язки рівняння є = f(хо). Якщо f(x) не дорівнює 0, за певним алгоритмом формується послідовність уточнень хn з використанням інформації про знак нев'язки, про значення самої нев'язки є або про швидкість її зміни де/дх.
   Вибір початкового значення х0 є важливим етапом, який істотно впливає на ефективність усієї процедури розв'язання і навіть на можливість одержання розв'язку.
   Знайомство з існуючими методами почнемо з прикладу простого одновимірного нелінійного рівняння. Припустимо, що це рівняння має не менше одного дійсного кореня.

4.1. Метод дихотомії
   Метод дихотомії (або метод поділу відрізка навпіл) передбачає послідовне обчислення значень функції в ряді точок. Перед використанням методу необхідно визначити відрізок, який містить лише один корінь рівняння. Для пошуку такого відрізка досить визначити точки a0 і b0, в яких знаки функції будуть протилежними, тобто f(а0)f(b0)<0, наприклад: f(а0)<0 і f(b0) > 0. Знайти такий відрізок можна, скориставшись графічним або табличним методом. Далі відрізок починають зменшувати, визначаючи на кожному кроці алгоритму координати його нових граничних точок an і bn за значеннями an-1, bn-1 та координату середньої точки попереднього відрізку:


У залежності від знаку функції в точці іщ новий відрізок [an, bn] функції встановлюється згідно з таким правилом:
     (4.2)

де n = 1, 2,..., mn - середня точка відрізку [аn-1, bn-1]

Рис. 4.1 Послідовне стискання інтервалу ізоляції кореня

Довжина інтервалу ізоляції кореня після виконання n кроків зменшується до
(bn-an) = 2-n(b0-a0)

а значення кореня а, обумовлене координатою середньої точки, і його похибки задаються виразами
     (4.3)

Із формули (4.3) випливає, що збіжність процесу обчислень дуже повільна, оскільки точність в одному десятковому розряді досягається за 3-4 кроки через те, що 10-1 приблизно дорівнює 2-3,3. Цей метод має абсолютну збіжність, тому ніяких вимог до виду і властивостей функції f(х) не існує (вони висуваються в разі використання інших більш швидкодіючих методів).

4.2. Метод простої ітерації
Метод простої ітерації полягає в тому, що рівняння (4.1) попередньо приводиться до канонічного вигляду:

     (4.4)

і ітерації виконуються за правилом
     (4.5)

причому задається початкове наближення х0.
Якщо процес обчислень збігається до розв'язку а рівняння (4.1), тобто а = ф(а), то в разі припущення, що функція ф(х) визначена і неперервна на відрізку, де треба знайти корінь, можна встановити взаємозв'язок між похибками обчислень на двох сусідніх ітераціях:
     (4.6)

Із виразу (4.6) випливає, що достатньою умовою збіжності методу простої ітерації, за якої |єn+і| буде менше |єn|, є умова
     (4.7)

Умова збіжіїості (4.7) ітерацій під час розв'язування нелінійних рівнянь є, власне кажучи, узагальненням достатньої умови збіжності ітераційного процесу під час розв'язання систем лінійних рівнянь, отриманої в розділі З у вигляді М < 1, якщо покласти |ф'*)| = М. Чим менше значення |ф'(x*)|, тим швидше збігається ітераційний процес.
У випадку, коли ф'(x*>0, похибки єn+і і єn будуть мати однакові знаки і збіжність хn+1 до а буде монотонною. Якщо ж ф'(x*) < 0, похибки єn+1 і єn матимуть різні знаки і наближення хn+1 буде збігатися до а, коливЮЧИСЬ біля а. Коли ф'(x*)>1, похибка єn+1 за абсолютним значенням більше єn, і наближення хn+1 буде відстояти від розв'язку а далі, ніж хn. Розв'язок а мов «відштовхує» наближення хn, близькі до нього, тому не буде збіжності послідовності хn до а.

Рис. 4.2 Можливі варіанти збіжності ітерацій: а - монотонна збіжність; б - коливальна збіжність; в - монотонна розбіжність; - коливальна розбіжність

Особливим випадком є умова, коли ф'(x*) = 0, і збіжність ітерацій замість лінійної, обумовленої виразом (4.6), стає нелінійною, зокрема квадратичною, як це буде показано в наступному розділі. Оцінювання глобальної похибки |хn+1 - а| зручно виконувати на основі значень, локальної похибки, тобто за значеннями наближень. Для цього додамо до правої частини виразу (4.6) і віднімемо від неї хn+1.

Після зведення подібних членів одержимо:

де |ф'(x*)| = M.
Якщо обчислювати похибку від початкового значення х0, то для поточної похибки на п-й ітерації згідно з виразом (4.6) справедливе співвідношення
     (4.8)

4.3. Метод Ньютона
Для прискорення збіжності ітераційного процесу чергове наближення xn+1 може бути визначене за формулою, де враховано як значення самої функції f(xn), так і швидкість її зміни f'(xn):

     (4.9)


Рис. 4.3 Ітераційний процес методу Ньютона

Цю формулу можна отримати безпосередньо з ряду Тейлора, в якому зберігаються два перших члени розкладу функції:
f(xn+1) = f(xn) + (xn+1 - xn)f'(xn) = 0,

звідки

де

Ефективність процедури (4.9) залежить від початкового наближення x0. наприклад, воно може збігатися з одним із кінців інтервалу ізоляції кореня [a,b], для якого виконується умова f(x)*f''(x) > 0.
Приклад 4.1


Знайдемо методом Ньютона розв'язок рівняння

для цього використаємо програму:

Вона дозволяє не тільки обчислити потрібну похідну f'(x) = 2sin(x + pi/3)cos(x + pi/3) - (x/2), але й отримати такі результати, якщо x0 = -1:

і якщо x0 = 2, де останній стовпець містить значення локальної похибки обчислень:

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

У припущенні, що


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

Поділивши співвідношення (4.10) на f'(x), з урахуванням формули (4.9) отримаємо

оскільки

одержуємо квадратичну залежність похибки на наступних ітераціях:
     (4.11)

Таким чином, метод Ньютона має квадратичну збіжність:
     (4.12)

Якщо з урахування вигляду функції вибрати початкове значення x0, так щоб на інтервалі

виконувалася умова

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

На практиці ж вибір значення x0 набагато складніший. Умови збіжності методу Ньютона, досліджені різними авторами, можна визначити у такий спосіб:
  • функція f(x) має бути принаймні двічі диференційована на інтервалі [a,b];
  • перша похідна функції не повинна перетворюватися на нуль f'(x) не дорівнює 0;
  • знаки функції на кінцях інтервалу мають бути різними sign f(a) не дорівнює sign f(b);
  • друга похідна на інтервалі [a,b] повинна зберегти свій знак f''(x) >> 0 (опукла) чи f''(x) < 0 (увігнута);
  • зростання функції (її крутість) має бути обмежене значенням

    що ускладнює роботу з експонентними функціями.

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

    Рис. 4.4 Приклади розбіжності методу Ньютона: а - функція з перегином; б - "зациклення" обчислень; в - експонентна функція

    Збіжність процедури обчислення кореня рівняння за методом Ньютона можна підвищити, якщо замість виразу (4.9) виконувати ітерації методу за формулою
         (4.13)

    і на кожному кроці обирати коефіцієнт демпфірування an з умови зменшення норми функції:
    ||f(xn+1)|| < ||f(xn)||.     (4.14)

    У формулі (4.13) співвідношення норм ||f(xn+1)|| і ||f(xn)|| залежить від значення коефіцієнта демпфірування an, і таку залежність можна вважати квадратичною:

    Рис. 4.5 Криві залежності K(an): 1 - розв'язання лінійної задачі за одну ітерацію amin = 1, K(1) = 0;

    У цьому можна пересвідчитися, проаналізувавши експериментальні криві Бреніна, ща наведені вище на зображенні, - незалежно від значення коефіцієнта відношення норм відхилу K(an) всі криві мають спільну дотичну в точці [0,1] із нахилом, рівним -2, і частки з оптимальними значеннями anmin.
    Спочатку на кожній ітерації оцінюється коефіцієнт K(1) за повного кроку збільшення

    і через точки [0,1] і [1, K(1)] проводиться парабола, коефіцієнти якої відповідно дорівнюють:

    Для параболи

    з урахуванням наведених значень коефіцієнтів обчислюються координати мінімальної точки amin і K(anmin):

    при цьому

    Крива 2 залежності коефіцієнта відношення норм відхилу K(an) відповідає вибору 0 < anmin < 1. Якщо ж одержуємо оцінку K(1) > 100 (крива 3), то параболічна апроксимація не виконується і коефіцієнт послідовно зменшується an = 1/2, 1/4, 1/8, ..., поки не виконається умова K(an) < 1. Якщо K(1) < 0,1, умова збіжності виконується з запасом, і обчислення проводяться з повним кроком збільшення аргументу an = 1.
    Процедуру послідовного підбору коефіцієнта демпфірування an без параболічної апроксимації для розглянутого вище випадку K(1) > 100 можна поширити і на всі інші випадки, коли в процесі обчислень не виконується умова (4.14), тобто ||f(xn+1)|| > ||f(xn)||. У цьому разі обирається деякий постійний коефіцієнт a0, наприклад a0 = 0, 75, що послідовно використовується на всіх ітераціях, для яких ||f(xn+1)|| > ||f(xn)||, при цьому поточне значення коефіцієнта дорівнює an = a0k, де k - кількість попередніх кроків, для яких не виконувалася умова (4.14).
    Починаючи з ітерації, для якої ця умова вже виконується, наявне значення коефіцієнта демпфірування an збільшується в 1/a0 разів на кожній ітерації аж до відновлення обчислень із повним кроком
    Приклад 4.2


    Розв'яжемо методом Ньютона рівняння f(x) = 10arctg(5 - x) - 1 = 0. У разі використання стандартної процедури методу (4.9) навіть за ненульового початкового значення ітерації не збігаються:

    Якщо перейти до процедури

    із вибором коефіцієнта an = 0,75 на кожному кроці обчислення за формулою an = a0k, то ітерації збігаються досить швидко:

    4.4. Метод Ньютона для кратних коренів
    Швидкість збіжності методу Ньютона падає, якщо рівняння f(x) має кратні корені. Метод Ньютона в цьому випадку втрачає квадратичну збіжність, і його збіжність стає лінійною. Швидкість збіжності залежить від кратності кореня q:

         (4.15)

    Наприклад, для q=2 збіжність методу Ньютона буде такою ж, як і збіжність методу дихотомії, тобто
         (4.16)

    Однак, якщо кратність кореня q відома заздалегідь, можна зберегти квадратичну збіжність методу Ньютона, модифікувавши основну формулу ітерацій:
         (4.17)

    У більшості випадків кратність коренів невідома, тому для збереження квадратичної збіжності на базі заданого рівняння з кратним коренем а формують інше рівняння
         (4.18)

    яке має одиничний корінь а незалежно від значення його кратності q у вихідному рівнянні f(x)=0. Це можна легко довести, якщо врахувати, що кратний корінь задовольняє як саме рівняння, так і похідні від нього, тобто f(j)(a)=0, якщо j

    переконуємося, що коли

    Тоді замість виразу (4.9) формула методу Ньютона для кратних коренів набуває такого вигляду:

    чи остаточно:
         (4.19)

    Приклад 4.3


    Знайдемо методом Ньютона кратні корені рівняння
    f(x) = x3 - 5x2 + 7x - 3 = (x - 3)(x - 1)(x - 1).

    Спочатку застосуємо стандартну процедуру:

    У разі вибору x0 = 0,5 корінь x1 = 1 локалізується з заданою точністю після шостої ітерації:

    Із застосуванням модифікованої процедури (4.19) за того ж початкового значення x0 = 0,5:

    невідомий корінь x1 = 1 локалізується після другої ітерації

    4.5. Обчислення коренів поліномів методом Ньютона
    Корені поліномів типу

         (4.20)

    можна знайти за допомогою методу Ньютона
         (4.21)

    при цьому для обчислення значень першої й другої похідних від поліноміальних функцій зручно використовувати рекурсивну процедуру Горнера для обчислення поліноміальних коефіцієнтів:
         (4.22)

    Знайдений корінь виключається з полінома, і порядок останнього зменшується:
         (4.23)

    Описана процедура повторюється n разів, поки не будуть виключені все корені. Однак часто поліноми мають комплексно-спряжені корені. У цьому випадку початкове значення вибирається також комплексно-спряженим zn=xn+iyn, і після знаходження пари таких коренів вони виключаються з полінома одночасно:

    Схема Горнера (4.22) при цьому змінюється:
         (4.24)

    Тут

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

    Приклад 4.4


    Знайдемо комплексно-спряжені корені полінома x2+1=0 для початкового комплексного значення x0=1+i. Скориставшись формулою методу Ньютона (4.9) проведемо обчислення:

    Вже після п'ятої ітерації переконуємось, що а = і.

    4.6. Визначення екстремальних точок функції методом Ньютона
    Задачу обчислення значень аргументу функції xe, за яких фуенкція f(xe) досягає своїх екстремальних значень (максимального чи мінімального), можна звести до задачі розв'язання нелінійних рівнянь, оскільки в даних точках похідна від функції дорівнює нулю, тобто f'(xe)=0. При цьому формула (4.9) набуває такого вигляду:

         (4.25)

    4.7. Точність обчислень у разі застосування методу Ньютона
    оскільки метод Ньютона є базовим у програмних комплексах математичного моделювання об'єктів і систем, розглянемо похибку обчислень даного методу. Для цього скористаємося відповідними співвідношеннями типу (4.8), виведеними раніше для методів простої ітерації (нелінійного й лінійного).
    Зведемо формулу (4.9) методу Ньютона до виду, характерного для методу простої ітерації:

         (4.26)

    де

    - похибка обчислень.
    Тоді зв'язок глобальної й локальної похибок визначається співвідношенням
         (4.27)

    де

    4.8. Метод Січних
    Застосовуючи формулу (4.9) чи (4.19) методу Ньютона, необхідно точно обчислювати похідні від розглянутих функцій, бажано аналітично. На практиці під час обчислень похідні апроксимують розділеними першими різницями:


    Використовуючи апроксимації у формулі методу Ньютона, можна перейти до методу січних:
         (4.28)

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

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

    Якщо ми маємо справу з кратними коренями, за аналогією з (4.19) можна записати:

    де u(x) = f(x)/f'(x).

    4.9. Метод хорд
    Об'єднання процедур методів січних і дихотомії дозволило побудувати метод хорд (або метод помилкової точки), в якому передбачено вибір початкових значень для ітерації xn і xn-1 з умови f(xn)f(xn-1)<0. При цьому кількість ітерацій у порівнянні з методом дихотомії зменшується за рахунок розбиття інтервалу ізоляції кореня не навпіл, а у відношенні f(a)/f(b). Чергове наближення визначається за формулою:

         (4.29)

    де с - нерухома точка.

    Рис. 4.7 Ітераційний процес методу помилкової точки

    За значення с береться той із кінців інтервалу локалызації [a,b], для якого f(x)f''(x)>0. Інший кінець інтервалу приймається за початкове наближення x0. Ітераційний процес закінчується в разі виконання умови

    де

    тому що існує оцінка

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

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

    f(x) = f(xn) + f(xn-1, xn)(x-xn) + f(x(sub>n-2, xn-1, xn)(x - xn-1)(x - xn),     (4.30)

    де застосовані перша і друга розділені різниці


    Рис. 4.8 Ітераційний процес методу Мюллера

    Ідея методу полягає в тому, що процес пошуку коренів рівняння (4.1) заміняється обчисленням коренів інтерполяційного многочлена (4.30), побудованого для f(x).
    Розв'яжемо рівння f(x) = 0. Позначимо різницю наближень на двох сусідніх ітераціях як z = xn+1 - xn і запишемо рівняння (4.30) для випадку x = xn+1 у такому вигляді:
         (4.31)

    де
         (4.32)

    Якщо продиференціювати рівняння (4.30):
    f'(x) = f(xn-1,xn) + f(xn-2, xn-1, xn)[(x - xn-1) + (x - xn)]

    і обчислити похідну для x = xn, можна переконатися, що
         (4.33)

    Розв'язуючи рівняння (4.31) отримаємо два кореня z(1) і z(2), на основі яких обчислюються корені x(1) і x(2). Як наступне наближення в методі Мюллера вибирається те значення x(1) і x(2), яке ближче до xn. Зокрема, для рівняння (4.31) знаходимо:

    Після очевидних спрощень із урахуванням (4.33) остаточно отримуємо:
         (4.34)

    Метод Мюллера зручний тим, що в загальному випадку дозволяє отримати комплексно-спряжені рівняння (4.1), використовуючи дійсні початкові наближення xn, xn-1, xn-2.