Умножим исходную матрицу A

92
Система (6.6) содержит n
2
уравнений и является линейной относительно искомых неизвестных x
ij
. Нетрудно заметить, что совокупность уравнений (6.6) состоит из n линейных систем, которые имеют одну и ту же матрицу коэффициентов при неизвестных (4.2).
Друг от друга системы уравнений отличаются векторами свободных членов
δ
ij
, причем в каждой системе только один свободный член равен единице (остальные равны нулю). Каждой j-й набор свободных членов определяет j-й столбец x
ij
(i = 1, …, n) искомой обратной матрицы A
−−−−1
Поиск элементов x
ij
обратной матрицы A
−−−−1
решением системы уравнений (6.6) эффективно осуществить алгоритмом Гаусса, который ранее успешно применялся в главах 4 и 5. Наличие нескольких столбцов свободных членов при неизменной матрице коэффициентов при неизвестных в системе уравнений (6.6) обусловливает некоторую специфику применения метода Гаусса к решению данной задачи обращения матрицы.
Прежде всего, строится расширенная матрица размером n
× 2n. К исходной матрице А приписывается справа еще единичная матрица Е того же порядка.
Далее к расширенной матрице (n – 1) раз применяется рекуррентное преобразование, аналогичное прямому ходу метода Гаусса в ходе решения системы линейных уравнений (см. § 4.3 и соотношения (4.17)).
В результате левая половина расширенной матрицы приводится к диагональному виду.
Затем реализуется процедура обратного хода метода Гаусса, причем на каждом шаге вычисления проводятся для всех столбцов свободных членов, что дает очередную строку искомой обратной матрицы.
В результате этих действий левая половина расширенной матрицы превратится в единичную, а в правой половине сформируется матрица
A
−−−−1
размером n
× n, обратная исходной А. В справедливости того, что правая половина преобразованной матрицы является обратной данной, можно убедиться непосредственным перемножением. Произведение полученной квадратной матрицы и исходной A дает единичную матрицу Е размером n
× n.
Пример 1. Используя метод Гаусса, найдем матрицу, обратную данной:

93
А
=












−
−
−
−
−
3 1
1 1
3 1
0 3
1 1
2 1
2 1
1 2

(6.7)
Сначала составим расширенную матрицу вышеописанного типа и поместим ее в таблицу 6.1. В крайнем левом столбце запишем порядковые номера строк. В следующие 4 столбца поместим исходную матрицу, в четырех правых столбцах – единичную матрицу размером 4
× 4. Элементы расширенной матрицы обозначим c
ij
, где i – номер строки, j – номер столбца.
Верхний левый элемент выбираем в качестве ведущего. Выделим его в таблице прямоугольником. Под штриховой линей в таблице 6.1 запишем строку, получаемую делением первой строки исходной расширенной матрицы на ведущий элемент. Присвоим полученной строке номер k=1. Очевидно, что первый элемент этой строки равен единице.
Таблица 6.1
i
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 1
2 3
4 2
1 3
1
−1 2
0
−1 1
−1
−1 1
2 1
−3 3
1 0
0 0
0 1
0 0
0 0
1 0
0 0
0 1
k
=1 1
−1/2 1/2 1
1/2 0
0 0
Затем строки с номерами i = 2, 3 и 4 исходной расширенной матрицы подвергаем преобразованию Гаусса, т.е. умножаем нижнюю строку с
номером k=1 на коэффициенты c
i
1
(i = 2, 3, 4) и вычитаем из соответствующей i-й строки. Иначе говоря, вычисляем элементы новой матрицы по формулам (4.17). Получаем новую матрицу из трех строк, у которой первый столбец с j=1 состоит из нулей. Запишем полученную матрицу в таблицу 6.2.
Теперь ведущим берем элемент c
22
= 5/2. Его также выделяем прямоугольником и делим на него строку с номером i = 2 в таблице 6.2.
Результат деления записываем ниже штриховой линии в таблице 6.2 в

94 строке, обозначенной k=2. Затем умножим полученную строку с номером k = 2 на коэффициенты c
32
и c
42
таблицы 6.2 и вычтем из соответствующих двух строк с номерами i = 3 и i = 4 той же таблицы.
Таблица 6.2
i
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 2
3 4
0 0
0 5/2 3/2
−1/2
−3/2
−5/2 1/2 0
−6 2
−1/2
−3/2
−1/2 1
0 0
0 1
0 0
0 1
k
=2 0
1
−3/5 0
−1/5 2/5 0
0
Получается матрица из двух строк, у которой уже два левых столбца состоят из нулей. Эту новую матрицу запишем в таблицу 6.3.
Таблица 6.3
i
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 3
4 0
0 0
0
−8/5 1/5
−6 2
−6/5
−3/5
−3/5 1/5 1
0 0
1
k
=3 0
0 1
15/4 3/4 3/8
−5/8 0
На следующем шаге ведущим элементом является c
33
=
−8/5. На него делится строка с номером i = 3 в таблице 6.3, а результат записываем под штриховой линией в таблице 6.3 в строке под номером k=3. Далее полученная строка умножается на коэффициент c
43
= 1/5 и вычитается из строки с номером i = 4 таблицы 6.3.
Результат вычитания
− единственная строка − помещается под номером i = 4 в таблицу 6.4. Наконец поделим строку c i = 4 таблицы 6.4 на коэффициент c
44
= 5/4 и результат запишем в ту же таблицу ниже штриховой линией под номером k=4.

95
Таблица 6.4

i
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 4
0 0
0 5/4
−3/4 1/8 1/8 1
k
=4 0
0 0
1
−3/5 1/10 1/10 4/5
Таким образом заканчивается прямой ход обращения матрицы.
Строки, пронумерованные индексом k, образовали расширенную матрицу, у которой левая половина приведена к диагональному виду
(см. табл. 6.5).
Таблица 6.5
k
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 1
2 3
4 1
0 0
0
−1/2 1
0 0
1/2
−3/5 1
0 1
0 15/4 1
1/2
−1/5 3/4
−3/5 0
2/5 3/8 1/10 0
0
−5/8 1/10 0
0 0
4/5
Таблица 6.5 содержит коэффициенты и свободные члены системы
(6.6), полученные преобразованиями прямого хода метода Гаусса. Числа в столбцах j = 1, …, 4 таблицы 6.5 можно рассматривать как коэффициенты при неизвестных элементах обратной матрицы, а числа в столбцах j = 5, …, 8 представляют собой четыре столбца свободных членов. Из уравнения (4.18) следует, что числа в строке k=4 и столбцах j
= 5, …, 8 составляют последнюю строку искомой обратной матрицы.
Для вычисления остальных строк следует выполнить операции обратного хода метода Гаусса.
Строка k=3 в таблице 6.5 эквивалентна системе уравнений
x
3j
+ (15/4) x
4j
= d
3j+4
,
(6.8) где j = 1,…,4, d
3j+4
– элементы строки k=3 и столбцов j = 1,…,4 таблицы
6.5.
Уравнения (6.8) аналогичны уравнению (4.19). Предпоследняя строка обратной матрицы получится решением уравнений (6.8), которое удобно выполнить формально с помощью таблицы 6.5. Для этого

96 достаточно из строки k=3 вычесть строку k=4, умноженную на элемент
d
34
=15/4. Таблица 6.5 преобразуется в таблицу 6.6.
Таблица 6.6
k
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 1
2 3
4 1
0 0
0
−1/2 1
0 0
1/2
−3/5 1
0 1
0 0
1 1/2
−1/5 3
−3/5 0
2/5 0
1/10 0
0
−1 1/10 0
0
−3 4/5
Элементы таблицы 6.6 в строке k=3 с j = 5,…,8 являются предпоследней строкой искомой обратной матрицы.
Следующий шаг обратного хода состоит в вычитании из строки k=2 таблицы 6.6 строки k=3, умноженной на элемент d
23
=
−3/5. При этом в правой половине строки k=2 сформируется вторая строка обратной матрицы А
−−−−1
. Результаты действий записаны в таблицу 6.7.
Таблица 6.7
k
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 1
2 3
4 1
0 0
0
−1/2 1
0 0
1/2 0
1 0
1 0
0 1
1/2 8/5 3
−3/5 0
2/5 0
1/10 0
−3/5
−1 1/10 0
−9/5
−3 4/5
На последнем шаге приходится решать уравнения
x
1j
+ d
12
x
2j
+ d
13
x
3j
+ d
14
x
4j
= d
1j+4
(j = 1,…,4),
(6.9) где d
kj
– элементы таблицы 6.7, а x
kj
(k = 2, 3, 4) – уже ранее вычисленные элементы строк 2, 3, 4 обратной матрицы. Для этого из первой строки таблицы 6.7 вычитается сумма остальных строк, каждый элемент которой предварительно умножается на d
1j
, где j – номер столбца (j=1,…,4).
В результате получается расширенная матрица, приведенная в таблице 6.8.

97
Таблица 6.8
k
j
=1
j
=2
j
=3
j
=4
j
=5
j
=6
j
=7
j
=8 1
2 3
4 1
0 0
0 0
1 0
0 0
0 1
0 0
0 0
1 2/5 8/5 3
−3/5 1/10 2/5 0
1/10 1/10
−3/5
−1 1/10
−1/5
−9/5
−3 4/5
Левая часть таблицы 6.8 содержит единичную матрицу, правая часть представляет собой искомую обратную матрицу А
−−−−1
:
А
−−−−1
=












−
−
−
−
−
−
5
/
4 10
/
1 10
/
1 5
/
3 3
1 0
3 5
/
9 5
/
3 5
/
2 5
/
8 5
/
1 10
/
1 10
/
1 5
/
2

(6.10)
Заметим, что перед расчетом А
−−−−1
не вычислялся детерминант исходной матрицы det A. Если бы det A = 0, то обратная матрица не существовала бы. В этом случае на одном из шагов прямого хода столбец, из которого берется очередной ведущий элемент, состоял бы только из нулей. Эта особенность указывалась в § 4.3 при описании метода Гаусса. Для независимой проверки вычислим детерминант матрицы (6.7) любым из способов, приведенных в предыдущей главе, и получим det A =
−10. Следовательно, матрица, обратная данной, существует.
Для проверки результата можно рассчитать обратную матрицу другим способом, описанным в начале параграфа, который выражает матрицу А
−−−−1
формулой (6.5). Однако вычислять 16 определителей 3-го порядка плюс еще один 4-го порядка – занятие трудоемкое.
Рациональнее провести перемножение полученной матрицы (6.10) и исходной (6.7). Произведение этих матриц оказывается равно единичной матрице, в чем может убедиться каждый читатель.
Следовательно, матрица (6.10) обратна данной (6.7), согласно определению (6.1).
В приведенном выше примере обращение матрицы демонстрировалось с помощью восьми таблиц. Столь подробное

98 изложение было проведено для увеличения наглядности используемого алгоритма Гаусса.
Использование компьютеров привело к распространению более эффективных методов. Например, алгоритм Гаусса
−Жордана позволяет провести обращение неособенной матрицы размером n
×n за n шагов.
Сначала строится расширенная матрица размером n
×2n, левая половина которой является исходной матрицей А, а правая половина
−−−− единичной матрицей.
На каждом k-м шаге алгоритма k-я строка расширенной матрицы делится на ведущий элемент. Элементы остальных строк (с номерами i
≠
k
) преобразуются по закону
)
(k
ij
c
=
)
1
(
−
k
ij
c
−
)
1
(
−
k
ik
c
)
(k
kj
c

(6.11)
Пример 2. Кратко проиллюстрируем алгоритм Гаусса-Жордана на матрице (6.7), обращенной в примере 1. Приведем только преобразуемые расширенные матрицы на каждом шаге алгоритма и соответствующие ведущие элементы:
B
0
=












−
−
−
−
−
1 0
0 0
3 1
1 1
0 1
0 0
3 1
0 3
0 0
1 0
1 1
2 1
0 0
0 1
2 1
1 2
,
c
11
= 2;
B
1
=












−
−
−
−
−
−
−
−
1 0
0 2
/
1 2
2
/
1 2
/
1 0
0 1
0 2
/
3 6
2
/
5 2
/
3 0
0 0
1 2
/
1 0
2
/
3 2
/
5 0
0 0
0 2
/
1 1
2
/
1 2
/
1 1
,
)
1
(
22
c
=5/2;

99
B
2
=












−
−
−
−
−
−
−
1 0
5
/
1 5
/
3 2
5
/
1 0
0 0
1 5
/
3 5
/
6 6
5
/
8 0
0 0
0 5
/
2 5
/
1 0
5
/
3 1
0 0
0 5
/
1 5
/
2 1
5
/
1 0
1
,
)
2
(
33
c
=
−8/5;
B
3
=












−
−
−
0 8
/
1 8
/
1 4
/
3 4
/
5 0
0 0
0 8
/
5 8
/
3 4
/
3 4
/
15 1
0 0
0 8
/
3 8
/
5 4
/
1 4
/
9 0
1 0
0 8
/
1 8
/
1 4
/
1 4
/
1 0
0 1
,
)
3
(
44
c
= 5/4;
B
4
=












−
−
−
−
−
−
5
/
4 10
/
1 10
/
1 5
/
3 1
0 0
0 3
1 0
3 0
1 0
0 5
/
9 5
/
3 5
/
2 5
/
8 0
0 1
0 5
/
1 10
/
1 10
/
1 5
/
2 0
0 0
1

Видно, что на 4-м шаге исходная расширенная матрица преобразовалась в матрицу, у которой левая половина является единичной матрицей, а правая половина
−−−− матрицей, обратной данной
(6.7), т.е. совпадающей с (6.10).
Для освоения других способов обращения матриц следует ознакомиться с более подробными курсами численных методов, указанных в списке литературы [1,3,16,17].

100
Глава 7. РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
§ 7.1. Выделение корней
В настоящей главе рассматриваются нелинейные уравнения, которые можно привести к стандартному виду:
f
(x) = 0,

(7.1) где f (x) – произвольная функция действительного переменного x.
Корнем уравнения (7.1) называется число
ξ, которое, будучи значением аргумента x =
ξ, обращает уравнение (7.1) в тождество.
Иначе говоря, при x =
ξ функция f (x) принимает нулевое значение:
f
(x =
ξ) ≡ 0.

(7.2)
Следовательно, проблема поиска корней сводится к отысканию нулей функции f (x). Выяснение наличия нулей функции f (x) и определение их количества целесообразно проводить, максимально используя свойства функции f (x) в левой части уравнения (7.1).
Например, если функция f (x) представляет собой полином n-й степени, то она не может иметь более n нулей.
Задача решения нелинейного уравнения вида (7.1) разбивается на две части. Первая часть заключается в нахождении интервалов, каждый из которых содержат внутри себя только один корень, т.е. интервалов, содержащих только один нуль функции f (x). Вторая часть состоит в вычислении значения корня с заданной погрешностью.
Первая часть является в принципе более сложной, так как не существует универсальных и эффективных алгоритмов выделения интервалов, содержащих единственный корень, для функции f (x) общего вида.
При поиске интервалов, содержащих единственный нуль функции
f
(x), находящейся в левой части уравнения (7.1), полезно опираться на следующие теоремы математического анализа.
Теорема 1. Пусть функция f (x) определена на интервале [a, b] и непрерывна на нем. Если на концах этого интервала функция f (x) имеет различные знаки, т.е.
f
(a)
⋅ f (b) < 0

(7.3)

101 то внутри интервала [a, b] находится хотя бы один корень уравнения
(7.1). Доказательство следует из непрерывности функции f (x) и тождества (7.2).
Замечание. В найденном интервале может находиться несколько корней: как нечетное, так и четное количество (см. рис.7.1).
-10
-5 5
10
x
-10
-5 5
y

Рис. 7.1. График знакопеременной функции.
На интервале ( –5, 10) данная функции имеет 3 нуля.
Теорема 2. Пусть функция f (x) определена на интервале [a, b] и непрерывна на нем. Если на концах этого интервала производная f
′(x) сохраняет знак и при этом для функции f (x) выполняется условие знакопеременности (7.3), то внутри интервала [a, b] находится единственный корень уравнения (7.1).
Доказательство от противного.
Таким образом, для нахождения интервалов, содержащих корни уравнения (7.1), целесообразно исследовать области знакопеременности функции f (x) и ее производной f
′(x).
Пример 1. Решить уравнение x
⋅ln(x) – 1 = 0.
Вычислим значения функции f (x) = x
⋅ln(x) – 1 в точках x =1 и x = e, где e – основание натуральных логарифмов. Получим следующие значения:
f
(1) = –1 < 0, f (e) = e – 1 > 0.

102
Производная f
′(x) = ln(x) + 1, т.е. положительна при x ≥ 1.
Следовательно, интервал [1, e] содержит один корень. Нетрудно доказать, что функция f (x) = x
⋅ ln(x) – 1 имеет единственный нуль в своей области определения.
При наличии современной компьютерной техники теоретический анализ иногда целесообразно сочетать с просмотром значений функции
f
(x) на интервале, границы которого задаются пользователем.
Отображение значений функции возможно как в текстовом, так и в графическом режиме. Если функция f (x) достаточно гладкая, то такой предварительный просмотр позволит довольно быстро выявить интервалы знакопеременности функции f (x), т.е. интервалы, содержащие корни уравнения (7.1). Для быстро осциллирующей функции или резко изменяющейся в некоторых диапазонах следует просматривать более узкие интервалы.
После того как найдены границы интервала, содержащего только один нуль функции f (x), следует применять какой-либо численный метод нахождения корня уравнения (7.1) на этом интервале. Наиболее простые и достаточно универсальные методы изложены в следующих параграфах настоящей главы. Эти методы работают достаточно быстро и устойчиво, но требуют определенного начального приближения x
0
искомого корня. Для получения начального приближения x
0
, как правило, необходимо знать границы интервала, содержащего искомый корень.
При постановке задачи решения нелинейного уравнения вида (7.1) задается допустимая погрешность
ε > 0 вычисления корня.
Следовательно, в процессе решения ищется не абсолютно точное, а приближенное значение корня. Приемлемым решением называется любое число x, которое удовлетворяет неравенству
x − ξ≤ε, где
ξ − точное значения корня заданного уравнения (7.1).
К сожалению, быстродействие современных компьютеров иногда провоцирует студентов на примитивный способ простого перебора. При этом пользователь вычисляет значения функции f (x) с шагом аргумента, равным величине погрешности
ε. Процессор современного персонального компьютера, совершающий миллионы операций в секунду, позволит даже таким возмутительно неэффективным способом найти интервал знакопеременности шириной
ε. В качестве искомого

103 корня берется середина этого интервала. Конечно, этот результат является приближенным значением корня с заданной точностью, но про подобные действия уже было сказано одним знаменитым человеком:
«Низкий класс, нечистая работа».

§ 7.2. Метод половинного деления
Простым, устойчивым и достаточно эффективным является метод
половинного деления
, называемый также бисекцией или дихотомией.
Пусть нам дано уравнение (7.1). Используя методику, изложенную в
§ 7.1, находим интервал знакопеременности [a, b], содержащий корень данного уравнения. Алгоритм вычислительного процесса бисекции изображен на рис. 7.2 в форме блок-схемы.
Принцип работы алгоритма очевиден из приведенной схемы, поэтому можно ограничиться краткими комментариями. В самом начале задается погрешность вычисления корня
ε и величина допустимой невязки
δ. Далее вводятся границы интервала [a, b], содержащего искомый корень, и проверяется свойство знакопеременности функции
f
(x) на концах заданного интервала. При отсутствии знакопеременности придется изменить границы начального интервала [a, b], в противном случае начинается итерационный процесс. На каждом шаге текущий интервал делится на две равные части точкой с. Эта точка является одной из границ нового, вдвое более узкого, интервала знакопеременности функции f(x). Процесс продолжается до тех пор, пока интервал знакопеременности не станет ýже заранее заданной погрешности вычисления корня
ε.
Кроме того, учитывается, что функция f(x) в области своего нуля может иметь очень большую производную. Тогда в найденной точке
c
≈ξ функция f(x) может иметь значение f(c), значительно отличающееся от нуля. Поэтому в алгоритм добавляется условие, из-за которого уменьшение интервала, содержащего корень, продолжается до тех пор, пока абсолютная величина
|f(c)| не станет меньше заранее заданного малого числа
δ.
За приближенное значение искомого корня принимается средняя точка последнего интервала знакопеременности.
Данный алгоритм легко программируется и успешно реализуется на современных компьютерах.

104










Рис. 7.2. Алгоритм схемы половинного деления.
a, b
– границы интервала, в котором происходит поиск корня,
ε и δ – заданные погрешности для корня и значения функции в точке найденного значения корня,
ξ – значение искомого корня.
Метод половинного деления дает последовательность приближений, сходящихся к точному значению корня, если функция f (x) непрерывна на предварительно выделенном интервале [a, b]. Следует заметить, что
ε
,
δ

с=(a+b)/2
a=c
b=c
a,b
(b>a)
f
(a)
.
f
(b)<0
f
(a)
.
f
(c)<0

(b
−
a <
ε)&(|f (c)|<δ)
ξ=c
да нет нет нет да да

105 монотонность функции f (x) на этом интервале, вообще говоря, не является обязательной.
Более быстрые алгоритмы, т.е. требующие меньшего количества вычислений для достижения заданной точности, рассмотрены в следующих параграфах этой главы.

§ 7.3. Метод Ньютона
Пусть на интервале [a, b] функция f (x) непрерывна, первая и вторая производные f
′(x) и f ′′(x) непрерывны и монотонны. Пусть x
n
– некоторое приближенное значение корня из интервала [a, b]. Тогда точное значение корня можно представить в виде
ξ = x
n
+ h
n
,

(7.4) где h
n
– малая поправка.
Разложим функцию f (x) в точке x = x
n
в ряд Тейлора по степеням малой величины h
n
до линейного члена:
f
(x
n
+ h
n
)
≈ f (x
n
) + h
n
⋅ f ′(x
n
).
(7.5)
Так как f (x
n
+ h
n
) = f (
ξ) = 0, то величина поправки выразится из уравнения (7.5) таким образом:
h
n
=
− f (x
n
) / f
′(x
n
).

(7.6)
Следовательно, очередное приближение корня можно выразить в виде рекуррентного соотношения:
x
n+
1
= x
n
− f (x
n
) / f
′(x
n
).

(7.7)
Исследования рекуррентного соотношения показали, что итерации по формуле (7.7) сходятся к истинному значению корня, если за начальное приближение
x
0
взять ту границу исходного знакопеременного интервала [a, b], для которой выполняется следующее условие:
f
(x
0
)
⋅ f ′′(x
0
) > 0

(7.8)

106
Вместо аналитического доказательства последнего утверждения приведем графические иллюстрации (см. рис. 7.3).





а





б
Рис. 7.3. Иллюстрация сходимости процесса Ньютона: а) f
′′(x) > 0 на всем интервале [a, b], б) f
′′(x) < 0 на всем интервале [a, b].

На рис. 7.3.а изображен случай, когда вторая производная f
′′(x) > 0 на всем интервале [a, b], f (a) < 0, f (b) > 0. Примем за начальное приближение корня x
0
= a. Величина производной f
′(x
0
) равна тангенсу
a
b
c
d
ξ
y
x
f
(x)
a
b
x
ξ
c
y
f
(x)

107 угла наклона касательной к графику функции f(x
0
) в точке x
0
Следовательно, первая поправка к корню (7.6) на рис. 7.3.а изобразится отрезком ad. Точка d (следующее приближение корня x
1
) лежит вне интервала [a, b]. Следовательно, дальнейшее использование итераций по формуле (7.7) некорректно.
Если за начальное приближение корня принять x
0
= b, то первый шаг процесса (7.7) даст следующее приближение корня x
1
= с, которое лежит внутри интервала [a, b]. При этом видно, что дальнейшие приближения
x
2
, x
3
, x
4
… будут располагаться слева от истинного значения корня
ξ.
Последовательность приближений x
n
сходится к
ξ справа налево вдоль числовой оси x.
На рис. 7.3.б приведен другой случай, когда на всем интервале [a, b] вторая производная f
′′(x) < 0. Аналогичные рассуждения вновь показывают справедливость условия сходимости (7.8). В данном случае последовательность приближений x
n
(n = 1, 2, …) сходится к истинному значению корня
ξ слева направо.
После выбора начального приближения корня запускается процесс
(7.7) и останавливается по достижении требуемой погрешности.
Ясно, что для использования метода Ньютона необходимо, чтобы производная f
′(x) нигде на интервале поиска корня не равнялась нулю.
Свойства рекуррентного соотношения таковы, что в общем случае условие
|x
n
− x
n
−1
|<ε не обеспечивает выполнения неравенства |x
n
−ξ|<ε.
В курсе численных методов получено следующее соотношение для оценки погрешности достигнутого приближения искомого корня x
n
:
| x
n
− ξ | ≤
1 2
2m
M
( x
n
− x
n
−1
)
2
,
(7.9) где M
2
– максимум абсолютной величины второй производной f
′′(x) на интервале [a, b]:
M
2
=
b
x
a
≤
≤
max
 f ′′(x),

(7.10)
m
1
– минимум абсолютной величины первой производной f
′(x) на том же интервале:
m
1
=
b
x
a
≤
≤
min
 f ′ (x).

(7.11)

108
Если, как и в предыдущем параграфе, задано максимально допустимое отклонение
δ значения функции от нуля в точке приближенного корня x
n
, то следует дополнительно проверять условие
| f (x
n
)
| < δ

(7.12) и при необходимости продолжать итерации по формуле (7.7).
В качестве искомого значения корня берется последнее приближение x
n
≈ ξ.
Важнейшим достоинством метода Ньютона является высокая скорость сходимости. Ясно, что количество шагов алгоритма Ньютона до достижения значения x
n+
1
≈ ξ существенно зависит от поведения производной f
′(x) на интервале [a, b]. Но практика показывает, что в большинстве случаев количество шагов алгоритма (7.7), по крайней мере, на порядок меньше, чем в методе половинного деления.
Недостатками метода
Ньютона являются следующие:
1) необходимость вычисления первой производной функции на каждом шаге процесса, 2) требование монотонности второй производной на интервале поиска корня, 3) необходимость вычисления второй производной на концах интервала [a, b].
Кроме того, применение данного метода становится «опасным» в тех случаях, когда вблизи корня первая производная близка к нулю. При этом очередная поправка к корню может быть большой и разложение
(7.5) до линейного члена станет неудовлетворительным. В этих случаях очередное значение корня может даже выйти из интервала [a, b], что приведет к расходимости процесса (7.7).
Примечание. Из рис. 7.3 видно, что последовательные приближения искомого корня получаются с помощью прямых линий, которые касательны к графику функция f (x). По этой причине метод Ньютона часто называется методом касательных.

§ 7.4. Метод секущих
Этот метод весьма сходен с методом Ньютона, но на каждом шаге итераций вместо производной используется разделенная разность
[f (x
n
)
− f (x
n
−1
)] / [x
n
− x
n
−1
].

(7.13)

109
Подстановка выражения (7.13) вместо производной в уравнение (7.7) дает рекуррентное соотношение для поиска корня уравнения (7.1) в следующем виде:
x
n
+1
= x
n
− f (x
n
)
⋅ [x
n
− x
n
−1
] / [f (x
n
)
− f (x
n
−1
)].
(7.14)
В качестве двух начальных приближений корня x
0
и x
1
удобно брать границы исходного интервала [a, b]. Иллюстрация применения метода секущих приведена на рис.7.4.









Рис. 7.4. Иллюстрация сходимости метода секущих.
ξ − точка истинного корня. В качестве начальных приближений взяты границы исходного интервала x
0
= a и x
1
= b. Следующим приближением корня является точка x
2
= c.
Для сходимости процесса (7.14) требуется лишь непрерывность функции f (x). Справедливости ради следует заметить, что при выделении интервала, содержащего единственный корень, ранее нам пришлось использовать монотонность первой производной f
′(x). Кроме того, ясно, что если на каком-либо шаге процесса (7.14) окажется f (x
n
) =
f
(x
n
−1
), то произойдет сбой алгоритма.
y
x
a
b
ξ
c
f
(x)

110
Условия сходимости, как правило, обеспечиваются на первом этапе решения уравнения (7.1) при выборе интервала знакопеременности функции f (x). При постоянстве знака второй производной f
′′(x) на интервале поиска корня процесс (7.14) сходится почти так же быстро, как и процесс Ньютона (7.7).
Для оценки погрешности достигнутого приближения корня можно воспользоваться очевидным неравенством:
| x
n
− ξ | ≤
1
)
(
m
x
f
n
,

(7.15) где величина m
1
определена формулой (7.11).
В курсе вычислительной математики доказано, что при использовании метода секущих выполнение условия
|x
n
− x
n
−1
|≤ε гарантирует неравенство
|x
n
− ξ|≤ε , т.е. достижение заданной погрешности вычисления искомого корня
ξ.
Скорость сходимости метода секущих, вообще говоря, несколько ниже, чем у метода Ньютона, что искупается гораздо более общими условиями сходимости. Это значит, что количество шагов процесса
(7.14) может быть несколько больше, чем у процесса (7.7). Но при использовании рекуррентной формулы (7.14) не приходится многократно вычислять значения первой производной f
′(x) на каждом шаге итерации и, можно вообще не рассматривать поведение второй производной f
′′(x) на интервале [a, b].
Пример 2. Для сопоставления методов Ньютона и секущих рассмотрим простейшее нелинейное уравнение x
2
− 2 = 0. Легко обнаруживается интервал, содержащий один корень: a=1, b=2.
Производная f
′(x)=2x монотонная на интервале [0, 1], вторая производная f
′′(x) является положительной константой. Процесс
Ньютона следует начинать с точки x
0
=1, в качестве начальных приближений метода секущих возьмем границы интервале [0, 1].
Результаты четырех первых итераций приведены в следующей таблице.

111
Таблица 7.1
Номер итерации n
x
n

(метод Ньютона)
x
n

(метод секущих)
1 1,500000000 1,333333333 2
1,416666667 1,400000000 3
1,414215683 1,414634146 4
1,414213452 1,414211438
Наблюдается быстрая сходимость обоих процессов к точному значению корня
ξ = 2 = 1,414213562… Метод Ньютона при том же числе шагов дает несколько меньшие отклонения от точного значения искомого корня
ξ по сравнению с методом секущих. С другой стороны, при использовании метода Ньютона приходится на каждом шаге итерации дополнительно вычислять значения производной f
′(x
n
).
Пример 3. Решение важной физической задачи о положении максимума спектра теплового излучения приводит к нелинейному уравнению вида:
5 e
−
X
+ X
– 5 = 0.

(7.16) где X – безразмерная переменная.
Корень X
0
уравнения (7.16) связан с физическими величинами следующим соотношением:
X
0
=
B
k
c
β
πh
2
,

(7.17) где h – постоянная Планка, k
B
– постоянная Больцмана, с – скорость света,
β – константа смещения Вина. Константа β связывает абсолютную температуру Т излучающего черного тела и длину волны
λ
max электромагнитного излучения, для которой спектр теплового излучения имеет максимум:
λ
max
T =
β.

(7.18)
Решив нелинейное уравнение (7.16), следует выразить из (7.17) параметр
β через найденный корень X
0
и вычислить значение размерной

112 константы смещения Вина. Эти операции предлагается выполнить читателям самостоятельно и сравнить результат с табличным значением, приведенным в любом добротном справочнике физических констант.

113
ГЛАВА

Скачать 0.96 Mb.

Поделитесь с Вашими друзьями: