69328

Методи розв’язування нелінійних рівнянь. Збіжність методів розв’язування нелінійних рівнянь

Лекция

Математика и математический анализ

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

Украинкский

2014-10-03

806 KB

5 чел.

Лекція 13. Методи розв’язування нелінійних рівнянь. Збіжність методів розв’язування нелінійних рівнянь

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

(6.1)

виконується переважно чисельними методами з властивими їм ітераціями і локальними апроксимаціями.

Розв’язок починається з вибору початкової точки x0 і оцінки нев’язки рівняння ε = f(x0) 0. Далі будується послідовність уточнень xi з використанням або інформації про знак нев’язки ( нев’язка змінює знак при проходженні через точку розв’язку), або про значення самої нев'язки ε, або про швидкість її зміни ∂ε/∂x.

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

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

6.1 Метод дихотомії (поділ відрізка навпіл)

Метод дихотомії (або метод поділу відрізка навпіл) припускає послідовне визначення значень функції і її знаків у ряді точок. Для існування відрізка, що містить хоча б один корінь рівняння (6.1), досить визначити точки a і b (рис. 6.1), для яких знаки нев’язки будуть протилежними : f(a)f(b)<0 .

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


 (6.2)

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

Після виконання n кроків інтервал ізоляції кореня зменшується до значення:

,

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

 (6.3)

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

Приклад 6.1.

Методом дихотомії розв’язати з точністю до 3% рівняння y = sin2(x + π/3) - (x/2)2. Підвищенню ефективності методу дихотомії на етапі початкового вибору інтервалу ізоляції кореня сприяє графічний метод, згідно з яким досить побудувати графік функції й оцінити приблизно значення коренів рівняння по абсцисах точок перетину графіка функції з віссю абсцис. Для нашого прикладу графік функції представлений на рис. 6.2.

Рис. 6.2. Графік функції з прикладу 6.1.

Вибираємо початковий інтервал ізоляції першого кореня I0 = [ - 1,0], тому що f( - 1) < 0, f(0) > 0

Результати послідовного стиску інтервалу представлені в табл. 6.1.

Таблиця 6.1. Результати застосування методу поділу відрізка I0 = [ - 1,0] навпіл

n

an

bn

mn

f (mn)

0

- 1

0

- 0.5

>0

1

- 1

- 0.5

- 0.75

<0

2

- 0,75

- 0,5

- 0,6250

>0

3

- 0,75

- 0,625

- 0,6875

>0

4

- 0,75

- 0,6875

- 0,7188

<0

5

- 0,7188

- 0,6875

- 0,7031

<0

Нев’язка після п'ятої ітерації дорівнює f (m5) = - 0,0098, а сам корінь локалізований з похибкою

.

Вибираємо початковий інтервал ізоляції другого кореня I0 = [1,2], тому що f(1) < 0, f(2) > 0 Результати послідовного зменшення інтервалу представлені в наступній табл. 6.2:

Таблиця 6.2. Результати застосування методу поділу відрізка I0 = [1,2] навпіл

n

an

bn

mn

f (mn)

0

1

2

1,5

<0

1

1

1,5

1,25

>0

2

1,25

1,5

1,375

<0

3

1,25

1,375

1,3125

>0

4

1,3125

1,375

1,3478

>0

5

1,3478

1,375

1,3594

<0

Нев’язка після п'ятої ітерації дорівнює f(m5) = - 0,0123.

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

6.2. Метод простої ітерації

Полягає в тому, що рівняння (6.1) попередньо приводиться до канонічного виду

 (6.4)

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

, (6.5)

причому задається початкове наближення x0.

Якщо процес обчислень збігається до розв’язку x* рівняння (6.1), тобто x* = φ(x*), то можна при припущенні, що функція φ(x)(x) визначена і неперервна на розглянутому відрізку, встановити взаємозв'язок між похибками обчислень на двох сусідніх ітераціях:

. (6.6)

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

. (6.7)

Умова збіжності (6.7) ітерацій при розв’язку нелінійних рівнянь є, власне кажучи, узагальненням достатньої умови збіжності (3.6) ітерацій при розв’язку лінійних рівнянь, отриманої в главі 3 у вигляді M < 1, якщо покласти |φ’(x*)| = M. Чим менше |φ’(x*)|, тим швидше збігається ітераційний процес.

При φ’(x*) > 0 похибки εn + 1 і εn будуть мати однакові знаки і збіжність xn + 1 до x* має монотонний характер ( рис. 6.3, а). Якщо ж φ’(x*) < 0 , похибки εn + 1 і εn мають різні знаки і наближення xn + 1 буде сходитися до x*, коливаючись біля x* (рис.6.3,б).

Коли |φ’(x*)| > 1, похибка εn + 1 за абсолютним значенням більше εn, і наближення xn + 1 буде відстояти від розв’язку x* далі, ніж xn Розв’язок x* мов «відштовхує» наближення xn, близькі до нього, у результаті чого не буде збіжності послідовності xn до x* (рис.6.3, в и г).

 

 

Рис. 6.3. Можливі варіанти збіжності ітерацій: а — монотонна збіжність, б — коливальна збіжність, в — монотонна розбіжність, г — коливальна розбіжність.

Випадок |φ’(x*)| = 0 є особливим випадком, коли збіжність ітерацій замість лінійної, обумовленим виразом (6.6), стає нелінійною, зокрема, квадратичною, як це буде показано в наступному параграфі.

Зручно при розрахунках оцінку глобальної похибки |xn + 1 - x*| вести за значеннями локальної похибки, тобто за значеннями, отриманих на сусідніх ітераціях за аналогією з формулою (3.8а).Для цього в праву частину виразу (6.6) додамо і віднімемо xn + 1:

.

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

,

де |φ’(x*)| = M.

Якщо врахування похибки вести від початкового значення x0, то для поточної nої ітерації згідно виразу (6.6) справедливе співвідношення:

. (6.8)

Приклад 6.3.

Нехай задана функція

f[x_] = x^(2.5] - Exp[ - x^(1.2)]

Для відділення дійсних коренів побудуємо графік у тій області, у якій, за нашими оцінками, повинні знаходиться потрібні нам корені. Для нашого прикладу x>0. При x = 0 значення функції дорівнює f(0) = - 1. З ростом x функція росте і, отже, знайдеться значення x*, при якому f(x*) = 0. Будуємо графік f(x)

Plot[f[t], {t,0,2}, AxesLabel - >{“x”, “f[x]”}]

Рис. 6.4.Графічний метод виділення кореня

Вихідне рівняння можна замінити наступним еквівалентним йому рівнянням, приведеним до виду (6.4):

.

Перевіримо для нього виконання умови (6.7). Визначимо в Mathematica функцію φ(x) = Exp( - 0,41,2) і знайдемо значення її похідної в околі кореня 

φ[_]: = E^( - 0.4^1.2]; φ1[x_] = Dt[φ[x],x]
Plot[φ1[x], {x,0,1}]

Рис. 6.5. Графік похідної рівняння φ(x) = Exp( - 0,41,2)

З рисунка випливає, що |φ’(x)| < 1, x[0,1] і, отже, процес ітерацій для x0[0,1] збігається. Оцінимо отриману похибку, використовуючи співвідношення (6.8), прийнявши в ньому q = M = max|φ’(x)|, x[0,1] і x1 = φ(x0):

, (6.9)

Якщо похибка обчислення кореня рівняння не повинна перевищувати наперед заданого значення ε, то згідно формулі (6.9) можна знайти необхідне число кроків обчислень n.

.

Знайдемо значення q для розглянутого прикладу

Q = - findMinimum[φ1[x], {x,0.5}]; q = q[[1]],
0.351148

Тепер можна знайти необхідне число ітерацій n, щоб забезпечити обчислення кореня з заданою точністю ε

x0 = 1.; ε = 0.0001; nε = Log[ε(1 - q)/Abs[φ[x0] - x0]]/Log[q]
8.15372

Округлимо отриманий результат до найближчого цілого зверху. Таким чином, обчислення по формулі (6.5) необхідно виконати дев'ять разів. При послідовному уточненні значення кореня рівняння та сама формула застосовується кілька разів для отриманих значень, у Mathematica це можна зробити, використовуючи оператор NestList[f,x,n], що дає перелік результатів застосування f до аргументу x, починаючи з нуля, загалом n разів. Виконаємо це в Mathematica

L = NestList[φ,x0,9]
{1., 0.67032, 0.78074, 0.742896, 0.75578, 0.751377,
0.752879, 0.752367 0.752542, 0.752482}

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

Ep = L[[10]] - L[[9]]
- 0.0000596766

Для збіжності велике значення має вибір функції φ(x). Помітимо, що не завжди вдається розв’язати рівняння (6.1) явно відносно x так, щоб привести його до виду (6.4) і задовольнити умові збіжності (6.7).

Приклад 6.4.

Нехай необхідно знайти розв’язок наступного рівняння

Для локалізації кореня побудуємо графік відповідної функції:

H[x_]: = xSin[x^2] - E^( - x/2); Plot [h[x], {x,0.5,1}]

Рис. 6.6. Графік функції xsin(x2) - e - x/2 = 0

Розв'яжемо задане рівняння відносно x, одержимо

.

Перевіримо для нього виконання умови (6.4). Визначимо в Mathematica функцію

і знайдемо значення її похідної в околі кореня, приблизне значення якого дорівнює 0.9:

ψ[x_]: = E^( - x/2)/Sin[x^2]; ψ1[x_] = Dt[ψ[x],x];
Plot[ψ1[x], {x,0.8,1}, AxesLabel - >{“x”,”ψ[x]”}]

Рис. 6.7. Графік значень похідної функції ψ(x) = e - x/2/sin(x2)

З приведеного графіка випливає, що в околі кореня |ψ(x)| > 1, тобто умова збіжності тут не виконана. Якщо виконати ітерації по формулі (6.5)

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

x0 = 0.85; L1 = NestList[ψ,x0,9]
{0.85, 0098867, 0.735723, 1.34347, 0.525148, 2.82473, 0.245467, 14.6883
0.000757322, 1.74291∙10
6}

У тих випадках, коли не вдається явно розв’язати вихідне рівняння (6.1) відносно невідомої x так, щоб в отриманому рівнянні (6.4) функція φ задовольняла умові збіжності (6.7), то функцію φ визначають у вигляді:

 (6.10)

причому допоміжна функція τ(x)не змінює знака на відрізку, де шукається корінь. Зокрема, якщо τ(x) = τ = const те одержимо метод релаксації

для якого φ(x) = 1 + τf(x) і метод сходиться за умови

 (6.11)

Якщо в деякому околі кореня виконуються умови

те метод релаксації збігається при τ(0,2/M1) Оптимальне значення параметра при цьому вибирається із співвідношення

. (6.12)

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

6.3 Метод Ньютона (дотичних)

Для прискорення збіжності ітераційного процесу чергове наближення xn + 1 може бути визначене по наступній формулі

. (6.13)

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

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

звідки ,

де

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

Приклад 6.5.

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

y = sin2(x + π/3) - (x/2)2 ,

розглянутого раніше в прикладі 6.1. Для цього обчислимо похідну

у' = 2sin(x + π/3)cos(x + π/3) - (x/2)

і заповнимо наступні табл.6.3 і 6.4, у яких останній стовпець містить значення локальної похибки обчислень.

Таблиця 6.3. Результати застосування методу Ньютона для знаходження першого кореня

n

xn

f (xn)

f ‘(xn)

Δxn

0

- 1

- 0,2478

0,5943

0,4109

1

- 0,5831

0,1154

1,0921

- 0,1057

2

- 0,6887

0,0045

1,0014

- 0,0045

3

- 0,6932

1,02*10 - 5

0,9969

- 1,02*10 - 5

Як видно з таблиці, похибка визначення кореня ε = 10 - 5досягнута за три ітерації.

Таблиця 6.4. Результати застосування методу Ньютона для знаходження другого кореня

n

xn

f (xn)

f ‘(xn)

Δxn

0

2

- 0,9911

- 1,1877

- 0,8345

1

1,1655

0,3019

- 1,5418

0,1958

2

1,3613

- 0,0155

- 1,6752

- 0,0093

3

1,3521

- 1,305*10 - 5

- 1,6723

- 7,8*10 - 5

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

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

 (6.14)

З співвідношення (6.14) поділом його на f(x) одержуємо з урахуванням формули (6.13):

Оскільки εn = a - xn і εn + 1 = a - xn + 1, то при xna одержуємо квадратичну залежність похибки на наступних ітераціях:

, (6.15)

чи квадратичну збіжність методу Ньютона:

. (6.16)

Якщо вибрати з урахуванням характеру функції (тобто значень f(a), f(b) початкове значення x0 так, щоб на інтервалі [a - ε0,a + ε0] виконувалася умова

,

то метод Ньютона завжди збігається, і збігається квадратично. Про швидкість збіжності можна судити по наступній табл.6.5 зміни похибки, побудованої для випадку mε00,9<1

Таблиця 6.5. Швидкість збіжності метода Ньютона

n

1

2

3

4

5

6

7

εn

0,81

0,66

0,44

0,19

0,036

0,0013

0,00016

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

1. функція f(x) повинна бути, принаймні, двічі диференційована на інтервалі [a,b];

2. перша похідна функції на інтервалі не повинна обертатися в нуль f(x)0;

3. знаки функції на кінцях інтервалу різні signf(a)signf(b);

4. друга похідна на інтервалі [a,b] зберігає свій знак: f(x) > 0 (опукла) чи f(x) < 0 (увігнута); 

5. зростання функції (її крутість) повинне бути обмежено значенням |f(x)/f(x)| (b - a), де
x = a при |f(a)| < |f(b)| і x = b при |f(a)| > |f(b)| , що ускладнює роботу з експонентними функціями, властивими моделям транзисторів.

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

 

Рис. 6.9. Приклади розбіжності метода Ньютона: а функція з перегином, б  випадок «зациклення» обчислень і в експонентна функція

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

 (6.17)

і вибирати на кожному кроці коефіцієнт «демпфірування» αn з умови зменшення норми нев’язки:

 (6.18)

Оскільки експериментально було показано, що на кожному кроці можна підібрати деяке оптимальне значення αn min, залежність частки норм ||f(xn + 1)|| і ||f(xn)|| від коефіцієнта «демпфірування» αn (рис. 6.10) зручно апроксимувати поліномом другого ступеня:

. (6.19)

Рис. 6.10. Криві залежності K(αn): 1 випадок розв’язку лінійної задачі за одну ітерацію αmin = 1, K(1) = 0; 2 випадок, коли 0.1  K(1) 100; 3 — випадок K(1) > 100

Спочатку на кожній ітерації виконується оцінка частки норм нев’язки K(1) при повному кроці збільшення Δxn і через точки [0,1] і [1,K(1)] з урахуванням, що всі криві рис.6.10 мають спільну дотичну в точці [0,1] з нахилом, рівним - 2, проводиться парабола (6.19), коефіцієнти якої відповідно рівні:

c0 = 1 з умови K(0) = 1;
c1 = - 2 з умови ;
c2 = 1 + M(1) з умови .

Для параболи (6.19) з урахуванням приведених значень коефіцієнтів визначаються значення αn min і K(αn min) з умови екстремуму функції:

,

при цьому

і 

Крива 2 розв’язку (рис. 6.10) відповідає вибору 0 <αn min < 1. Якщо ж одержуємо оцінку K(1) > 100 (крива 3), то параболічна апроксимація не виконується і коефіцієнт демпфірування послідовно зменшується αn = 1/2,1/4,1/8,…, поки не виконається умова K(αn) < 1 При K(1) < 0,1 умова збіжності виконується з запасом і обчислення проводяться з повним кроком збільшення аргументу αn = 1.

Процедуру підбору коефіцієнта демпфірування αn без параболічної апроксимації, проведену вище для випадку K(1) > 100 (крива 3), можна, з метою спрощення обчислень, поширити і на всі випадки, коли в процесі обчислень не виконується умова (6.18), тобто ||f(xn + 1)|| > ||f(xn)||. З цією метою обирається деякий постійний коефіцієнт α0, наприклад, α0 = 0,75, що послідовно застосовується на всіх кроках обчислень, для яких ||f(xn + 1)|| > ||f(xn)||, при цьому

, (6.20)

де k — число попередніх кроків, на яких не виконувалася умова (6.18).

З моменту часу, коли ця умова починає виконуватися, наявне значення коефіцієнта демпфірування αn поступово збільшується в 1/α0 раз на кожному кроці аж до відновлення повного кроку.

Приклад 6.6.

Розв’язати методом Ньютона рівняння

(x) = 10*ATG(5 - x) - x = 0

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

Таблиця 6.6. Результати обчислень згідно стандартноюпроцедури (6.13)

x

y

0

13,734008

9,9190055

- 23,621369

- 6,9911018

21,867039

13,463214

- 27,995047

- 11,143668

26,232984

14,123512

- 28,739764

- 11,566567

26,671636

14,170702

- 28,792527

- 11,594276

26,70035

14,173701

- 28,795879

- 11,596027

26,702164

14,17389

- 28,79609

Якщо перейти до процедури (6.17) з вибором коефіцієнта «демпфірування» αnна кожному кроці обчислення по формулі (6.20) , то ітерації збігаються і досить швидко (табл.6.7).

Таблиця 6.7. Результати обчислень згідно модіфікованої процедури (6.20)

x

y

αn

0

13,734008

1

9,9190055

- 23,621369

1

- 2,7635749

17,190525

0,75

5,5493886

- 10,573125

0,5625

5,035591

- 5,3913511

0,421875

4,7595795

- 2,400154

0,5625

4,5873784

- 0,67398327

0,75

4,5167682

- 0,015335687

1

4,5150843

- 9,0085024E - 6

1

До речі, вихідне рівняння y(x) = 10∙atg(5 - x) - x = 0 при x0 = 0 пакет Mathematica розв’язати не може.

6.3.1. Метод Ньютона для кратних коренів

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

 (6.21)

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

. (6.22)

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

. (6.23)

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

 (6.24)

яке має одиночний корінь а незалежно від значення його кратності q у вихідному рівнянні f(х). Це можна легко довести, якщо врахувати , що кратний корінь задовольняє як саме рівняння, так і похідні від нього, тобтоf(j)(α) = 0 при j < q. Записавши для розглянутого випадку розклад в ряд Тейлора для функції і її похідної:

,

переконуємося, що при x→∞ lim[y(x)/(x - α)] = 1/q.

Тоді замість виразу (6.13) формула ітерацій Ньютона для кратних коренів приймає вид:

 

чи остаточно:

 (6.25)

Приклад 6.7.

Знайти, користуючись методом Ньютона, кратні корені рівняння

.

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

,

При виборі x0 = 0 корінь x1 = 1 локалізується після шостої ітерації (табл..6.8).

Таблиця 6.8. Результати обчислень згідно стандартною процедури (6.13)

i

xi

0

0

1

0,428571

2

0,685714

3

0,832865

4

0,913329

5

0,955783

6

0,977655

Тепер застосуємо модифіковану процедуру (6.25) при тому ж початковому значенні x0 = 0:

.

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

Таблиця 6.9. Результати обчислень згідно модіфікованої процедури (6.25)

i

xi

0

0

1

1,105263

2

1,003081

3

1,000003

Точність обчислень по методу Ньютона

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

Приведемо формулу (6.13) методу Ньютона до виду, який властивий методу простої ітерації:

, (6.32)

де

 

δ похибка обчислень.

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

 (6.33)

де, як і раніше у формулі (6.16),|kΔxn|m ≤ 1 і k ≥ |f(x)/f(x)|.

Приклад 6.9.

Оцінити похибку значень табл. 6.3 прикладу 6.5 для першого кореня рівняння

y = sin2(x + π/3) - (x/2)2,

у який x3 = - 0,6887 і Δx3 = - 0,0045 + 0,5∙10 - 4.

Перша і друга похідні від розглянутого рівняння рівні:

y’ = 2sin(x + π/3)cos(x + π/3) - (x/2);

Звичайно, можна було б за допомогою Mathematica оцінити максимально можливі значення цих похідних для визначення коефіцієнта k у виразі (6.33), як це б зроблено в прикладах 6.2 і 6.3. Але можна зробити це і наближено, приймаючи x = - 0,7 і вибираючи максимально можливі значення тригонометричних функцій, у результаті отримаємо: |f(x)| < 1, |f(x)| < 3/2.

Тоді

 

6.3.Метод січних (хорд)

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

.

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

(6.34)

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

 

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

Для випадку кратних коренів за аналогією з (6.25) можна записати

,

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

Приклад 6.7.

Розв'язати методом січних рівняння

y = sin2(x + π/3) - (x/2)2 ,

розглянуте вище в прикладах 6.4 і 6.5 і розв’язане методами дихотомії і методом Ньютона.

Для знаходження першого кореня обираємо два початкових значення x0 = - 1,1 і x1 = - 1.

Результати наступних ітерацій зведені в таблиці 6.11. Для локалізації другого кореня змінюємо початкові значення на x0 = 2 і x1 = 2,1, таблиця. 6.12.

Таблиця 6.11. Результати застосування методу січних для знаходження першого кореня

n

xn

f (xn)

Δxn

0

- 1,1

- 0,2997

1

- 1

- 0,2478

0,4772

2

- 0,5228

0,1824

- 0,2023

3

- 0,7251

- 0,0312

0,0296

4

- 0,6955

- 0,0023

0,0023

5

- 0,6932

Таблиця 6.12. Результати застосування методу січних для знаходження другого кореня

n

xn

f (xn)

Δxn

0

2

- 0,9911

1

2,1

- 1,1025

- 0,9899

2

1,1101

0,3857

0,2565

3

1,3666

- 0,0244

- 0,0153

4

1,3514

0,0012

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

Її можна оцінити, повторивши викладки, зроблені при виведенні формули (6.15). Отримаємо вираз

(6.35)

який можна апроксимувати формулою

 (6.36)

при цьому  (6.37)

Таким чином, збіжність методу хорд залишається нелінійною з показником 1,618 замість 2, як у методу Ньютона.

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

 (6.38)

і для якого

(6.39)

6.3.1. Метод хорд (або метод помилкової точки)

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

, (6.40)

де c – нерухома точка.

  

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

У якості c береться той з кінців вихідного інтервалу ізоляції [a,b], у якому f(x)f(x) > 0 .

Інший кінець відрізка приймається за початкове наближення x0. Ітераційний процес закінчується при виконанні умови

,

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

(6.41)

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

6.3.2.Комбінований метод

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

.

 

Рис. 6.13. Ітераційний процес комбінованого методу

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

, (6.41)

де xn(k) і xn(c) — наближення кореня, отримані відповідними методами.

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

Приклад 6.8.

Знайти корінь рівняння, заданого в прикладі 6.2, комбінованим методом. Знайдемо початкове наближення для методу дотичних, виходячи з виконання умови (6.13а). Нехай
x0 = 1. Визначимо першу і другу похідні функції

f1[x_] = Dt[f[x],x]; f2[x_] = Dt[f[x],{x,2}]; f2[1]*f[1]

2.0914

Отримане значення позитивне, отже, умова виконана.

Для методу хорд у якості початкового наближення необхідно задати два значення x0 і x1 таких, щоб була виконана умова f(x0)f(x1) < 0. Виберемо, наприклад, x0 = 1, x1 = 0,5, тоді остання умова буде виконано. Обчислення послідовних наближень виконується за формулами (6.13) і (6.16)

z1 = 0.5; z0 = 1; x = 1.; k = 0; ep = 0.0001; ek = 10;
While[ek > ep, z1 = z1-f[z1]*(z0-z1)/(f[z0]-f[z1]);
x = x-f[x]/f[x1]; ek = Abs[z1-x]; z0 = x; k = k+1];
Print[stringForm[“x = ’ ’, z1 = ’ ’, ek = ’ ’, k = ’ ’, x, z1,ek,k]] x = 0.752497325379611226, z1 = 0.752496715923983128, ek = 6.09455628075572519

Отримані результати показують, що за три ітерації знайдено значення шуканого кореня з похибкою

.

Розв’язок останнього рівняння можна отримати наступним оператором пакета Mathematica, що дає той же самий результат

FindRoot [f[u] = = 0, {u,1}]

6.4.Метод Мюллера

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

(6.42)

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

 

 

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

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

 (6.43)

де

. (6.44)

Якщо продиференціювати рівняння(6.42) :

і оцінити отриману похідну при x = xn , можна переконатися, що 

(6.45)

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

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

. (6.46)

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

Приклад 6.9.

Знайти корені наступного рівняння f(x) = Sin[x] + Cos[x] - 2 = 0.

Побудуємо графік функції y = f(x) 

 

Рис. 6.15. Графік функції з прикладу 6.9

Як видно з графіка, дійсних коренів розглянуте рівняння не має (графік не перетинає вісь абсцис). Виберемо три довільних значення для аргументу, наприклад, x = 0, 0,75, 1,5. Складемо таблицю значень функції в цих точках

T = Table[{x, f[x]}, {x, 0, 1.5, 0.75}]

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

P = InterpolatingPolynonial[T,x]

Знайдемо нулі отриманого багаточлена

Z = Solver [P =  = 0, x]

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

A = x/.Z[[1]]; b = x/.Z[[2]];

Внесемо один з них, наприклад a, у таблицю T. Спочатку видалимо перший елемент списку T, а потім допишемо рядок а.

T1 = Rest[T]; T = Append[T1, {a, f[a]}]

Повторюючи обчислення наближеного значення кореня, отримаємо

Проробивши те ж саме ще раз, отримаємо більш точне значення кореня. Знайдемо корені цього ж рівняння в Mathematica і порівняємо з отриманими вище значеннями 

Nsolver [f[x] =  = 0, x]

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

6.5. Методи розширення області розв’язку

6.5.1. Метод Давиденка продовження розв’язку по параметру

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

Відповідно до методу продовження розв’язку по параметру початкова, важко розв'язувана задача f(x) = 0, перетворюється в послідовність задач, що залежать від параметра, який вводиться в рівняння:

   (6.47)

параметр k вибирається так, щоб

де g(x) = 0 - легко розв'язувана задача.

У результаті розв’язок задачі є функцією від введеного параметра x = (k).

Включення параметра k у розв'язуване рівняння може бути виконане по-різному. Наприклад, можна скласти композицію функцій f(x) і g(x):

(6.48)

Початкове наближення на наступній ітерації також можна вибирати різним чином:

1 припасовуванням x0(ki) = x(kk - 1);

2 інтерполяцією  (6.49)

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

Є й інша можливість розв’язати рівняння (6.47), диференціюючи його по параметру k і приводячи його до еквівалентного диференційного рівняння типу

,

звідки

. (6.50)

Методи розв’язку звичайних диференціальних рівнянь, розглянуті в главах 8-10, не гарантують розв’язку рівняння (6.50) при зміні параметра k від 0 до 1 більш ефективно, порівняно з дискретною зміною значень k і використанням виразу (6.49).

При моделюванні електронних схем, наприклад, дискретний параметр k, що вводиться в їхні рівняння, дозволяє змінювати напругу загального джерела живлення від малих значень до номінального, якщо розв’язок задачі статики при номінальному значенні напруги не сходиться через наявність у рівняннях схеми експоненціальних рівнянь функцій p-n переходів транзисторів, як це зроблено в програмах PSPICE, ALLTED [20].

Приклад 6.10.

Розв’язати, використовуючи метод Давиденка при x0 = 0 експоненціальне рівняння

f(x) = x + (EXP(x/0,25) - 1) - 50 = 0.

Спосіб прямо розв’язати це рівняння методом Ньютона виявляється неуспішним, оскільки за припустиме число ітерацій збіжність обчислень не досягається (таблиця 6.13.):

Таблиця 6.13. Результати рішення експоненціального рівняння методом Ньютона

x

f(x)

0

- 50

10

2,3538527E17

9,75

8,65934E16

9,5

3,1855932E16

9,25

1,1719142E16

9

4,3112315E15

8,75

1,5860135E15

8,5

5,8346174E14

8,25

2,1464358E14

8

7,896296E13

7,75

2,904885E13

7,5

1,0686475E13

7,25

3,9313343E12

7

1,4462571E12

6,755,3204824E11

5,3204824E11

6,5

1,9572961E11

6,25

7,2004899E10

6

2,6489122E10

5,75

9,7448034E9

5,5

3,5849128E9

Тому вихідне рівняння приводиться до виду (6.47) введенням параметра k:

F(x,k) = f(x) = x + (EXP(x/0,25) - 1) - 50k - 5(1 - k) = 0.

виборі k = 0 розв’язок отримуваного рівняння методом Ньютона

F(x,0) = x + (EXP(x/0,25) - 1) - 5 = 0

сходиться за 7 ітерацій при x01 = 0 (таблиця 6.14.):

Таблиця 6.14. Результати рішення експоненціального рівняння при виборі k = 0

x

F(x,0)

0

- 5

1

49,59815

0,7739297

16,877051

0,58517473

4,9736697

0,46829566

0,97727493

0,43214838

0,064874396

0,42939139

0,00034125817

Тепер вибираємо k = 0,4 і при x02 = 0,42939139 розв’язуємо рівняння

F(x,0,4) = x + (EXP(x/0,25) - 1) - 23 = 0.

Обчислення сходяться також за 7 ітерацій (таблиця 6.15.):

Таблиця 6.15. Результати рішення експоненціального рівняння при виборі k = 0,4

x

F(x,0,4)

0,42939139

- 17,999659

1,2024465

99,907778

0,99930809

31,446559

0,85557802

7,4957536

0,79491348

0,83334593

0,78633587

0,013988626

0,78618692

4,1222073E - 6

І нарешті, при k = 1 повертаємося до вихідного рівняння

F(x,1) = f(x) = x + (EXP(x/0,25) - 1) - 50 = 0

і розв’язуємо його при x03 = 0,78618692. У результаті отримуємо шуканий розв’язок за 5 ітерацій (таблиця 6.16.):

Таблиця 6.16. Результати рішення експоненціального рівняння при виборі k = 1

x

f(x)

0,78618692

- 26,999996

1,0738638

23,439481

0,99426294

3,3537477

0,97862322

0,10227039

0,97811567

0,00010323141

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

Правильність отриманого розв’язку підтверджується результатом розв’язку вихідного рівняння за допомогою пакета Mathematica, у якому реалізований метод Ньютона з коефіцієнтом демпфіруючого типу (6.17) і (6.20):

In[8]: = FindRoot[P =  = 0, x]
Out[8] = 0.978115

 

Рис. 6.16. Графік функції f(x) = x + (EXP(x/0,25) - 1) - 50

6.5.2. Метод пошуку кривої розв’язку

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

 

Рис. 6.17. Ітерація в методі пошуку розв’язку

У початковій точці з координатою  будується дотична, як і в методі Ньютона, але на цю дотичну опускається перпендикуляр із точки, що лежить на осі x, і перетинає дотичну в точці С і криву функції в точці В. Координата точки В приймається за наступне наближення
xi + 1, у ній проводиться нова дотична і т.д. Для знаходження значення xi + 1 використовуються ітерації по стандартному методу Ньютона, початкове значення xi + 10 для якого задається координатою точки перетину перпендикуляра з дотичною. Можна ввести додатковий параметр λ, за допомогою якого контролюють кут перетинання додаткової прямої з дотичною, змінюючи його від прямого кута (λ0 = 1) до збігу цієї лінії з віссю xλ0 = 0, що відповідає звичайному методу Ньютона.

Координата точки С знаходиться з умови перетину двох перпендикулярних прямих:

або

,

звідки

або . (6.51)

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

,

по якому будується рівняння

;

розв'язуване методом Ньютона:

З урахуванням значення похідної

остаточно отримуємо:

(6.52)

Приклад 6.11.

Розв’язати рівняння f(x) = 3 arctg (3 - х) = 0 методом Ньютона і методом пошуку кривої розв’язку і порівняти отримані результати при x0 = 0.

Графік використовуваної функції показаний на рис. 6.18.

 

Рис. 6.18. Графік функції з прикладу 6.11

Організовуючи обчислення за формулою (6.13) методу Ньютона отримуємо результати зведені в таблицю 6.17. з таблиці видно, що розвязок не збігається:

Таблиця 6.17. Результати обчислень за методом Ньютона

x

f(x)

0

3,7471373

12,490458

- 4,3974441

- 120,99951

4,6881959

23908,94

- 4,7122635

- 8,9765283E8

4,712389

- 2,5164787E36

4,712389

9,947327E72

4,712389

Якщо задане рівняння розв’язувати, використовуючи формули (6.51) і (6.52) при λ = 1 і x0 = 0, то отримуємо розв’язок, який швидко сходиться за три ітерації, (таблиця 6.18.):Якісно характер розв’язку ілюструється на рис. 6.19.

Таблиця 6.18. Результати обчислень за методом пошуку кривої розв’язку

x

f(x)

0

3,7471373

2,505267

1,3782754

2,9931483

0,020554758

2,999315

0,0020551023

 

Рис. 6.19. Ітераційний процес методу пошуку кривої розв’язку

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

PAGE  124


 

А также другие работы, которые могут Вас заинтересовать

11473. Культура і етнос. Ментальне поле культури. Проблема культурної ідентичності 460.5 KB
  Лекція 7. Культура і етнос. Ментальне поле культури. Проблема культурної ідентичності. Термін етнос у давньогрецькій мові має більш як десять значень: народ племя юрба стадо стан соціальна група клас тощо. Етнічна спільність це спільнота повязана певною ...
11474. Роль культурних орієнтацій у розвитку суспільства 650.5 KB
  Лекція 8. Роль культурних орієнтацій у розвитку суспільства. Поняття культурна орієнтація трапляється все частіше в наукових і публіцистичних творах гуманітарного спрямування. Воно постійно використовується в повсякденному спілкуванні. Як і кожне багатозначне понят
11475. Глобалізація і культура 91 KB
  Лекція 9. Глобалізація і культура. Глобалізаційні процеси в сучасній світовій культурі та теорія модернізації Питання взаємовпливу культур різного типу як проблема теорії культури була поставлена ще культурною антропологією початку XX ст. Значний внесок у її роз
11476. Семіотика культури 59 KB
  Лекція 10. Семіотика культури Семіотика культури передбачає вивчення соціокультурних процесів і явищ як системи висловлювань культурних текстів культурних мов яке здійснюється насамперед лінгвістичними методами. Засновник семіотики американський фі...
11477. ТЕХНІКА, КУЛЬТУРА, ПРИРОДА ЛЮДИНИ 209 KB
  ЛЕКЦІЯ 11 ТЕХНІКА КУЛЬТУРА ПРИРОДА ЛЮДИНИ Питання про місце техніки в соціальній історії в житті сучасного суспільства в життєвих цінностях та орієнтирах людини є одним із найбільш болючих. Ціла низка гуманістично налаштованих мислителів вважає що техніка з ї...
11478. КУЛЬТУРА І ПОЛІТИКА 182.5 KB
  Лекція 13. КУЛЬТУРА І ПОЛІТИКА Політична культура являє собою динамічну та водночас відносно усталену систему політичних цінностей та орієнтацій моделей поведінки характерних для певної держави суспільства цивілізації. Від її стану харак
11479. ХУДОЖНЯ КУЛЬТУРА ПЕРШОЇ ПОЛОВИНИ XX ст.: АВАНГАРДНІ СТИЛІ ТА НАПРЯМИ 314.5 KB
  Лекція 14.ХУДОЖНЯ КУЛЬТУРА ПЕРШОЇ ПОЛОВИНИ XX ст.: АВАНГАРДНІ СТИЛІ ТА НАПРЯМИ Проблема сприйняття та розуміння художнього авангарду Новий зміст у традиційній мистецькій формі фовізм кубізм експресіонізм супрематизм Становлення нових форм мистецтва: футуризм ...
11480. ХУДОЖНЯ КУЛЬТУРА ДРУГОЇ ПОЛОВИНИ XX ст 159.5 KB
  Лекція 15. ХУДОЖНЯ КУЛЬТУРА ДРУГОЇ ПОЛОВИНИ XX ст. Термінологічна картина художнього авангарду Візуальні та інформаційні мистецтва: гіперреалізм концептуалізм бодіарт лендарт та ін. Акціонізм хепенінг перформанс рольова гра За відомим виразом ми живе
11481. ІДЕЯ РІВНОПРАВНОСТІ КУЛЬТУР У СУЧАСНОМУ СВІТІ 556 KB
  Лекція 16. ІДЕЯ РІВНОПРАВНОСТІ КУЛЬТУР У СУЧАСНОМУ СВІТІ Ідея про рівноправність культур формується на Заході починаючи з античної культури у безпосередньому звязку з розвитком поняття про Закон і його роль в суспільстві. Антична Греція як країна класичної демократі