Часть 45

Численные методы

Эта глава стоит особняком от всего остального, т.к. мы рассмотрим некоторые практические приложения математического анализа, а именно численные методы, но не претендуя при этом на строгость изложения.

1Метод бисекции

Начнем с самого простого метода. Пусть есть функция ( ): [ , ] →R, причем ( ) C[ , ]. Предположим, что ( ) 6 0, ( ) > 0 (или наоборот). Данный метод позволяет найти какое-либо одно решение уравнения ( ) = 0 на отрезке [ , ]. Не стоит путать этот метод с двоичным поиском! Двоичный поиск требует монотонности функции на отрезке [ , ] и используется, в основном, при поиске нулей у функций ( ): Z → R, т.е. не обладающих непрерывностью.

Теперь подробнее о методе. Для начала заметим, что теорема 13.4 гарантирует существование хотя бы одного решения. Найдем его следующим образом:

1.Проверим сначала, что ( ) = 0 или ( ) = 0. Если хотя бы одно из условий истинно, то искомый корень, очевидно, найден.

2.Примем = +2 . Теперь возможно 3 варианта:

( ) = 0. В этом случае точка — искомый корень.

( ) > 0. Тогда если ( ) < 0, то поставим = , иначе = .

( ) < 0. Тогда если ( ) < 0, то поставим = , иначе = .

3.Перейдем к шагу 2.

Проанализируем данный метод. Суть условий в пункте 2 заключается в том, чтобы сохранить инвариант ( ) 6 0, ( ) > 0 (или наоборот). Как несложно убедиться, данный инвариант действительно сохраняется. После каждой итерации алгоритма новые вычисленные значения

и находятся уже ближе к искомому корню. Заметим, что данный метод может никогда не сой-

√

тись. Например, если рассмотреть ( ) = 2− на отрезке [0, 2]. В этом случае каждое и является рациональным, в то время как искомое решение иррационально. Последовательность получаемых приближений составляет систему стягивающихся сегментов, причем

lim = lim = √2.
→∞	→∞

Однако численные методы призваны находить приближенные решения. Здесь решение очевидно, но, как мы увидим позднее, честное решение некоторых уравнений (или нахождение каких-либо других величин) крайне трудоемко. Чтобы искать приближенное решение, в нашем случае можно прибегнуть к использованию удовлетворяющей нас абсолютной погрешности. Иными словами, мы хотим найти такое , что | ( )| < . Соответственно, именно таким образом и следует производить проверку «равенства» ( ) = 0.

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

√

найдется такое , что | ( )| < и | ( )| < . Например, пусть необходимо вычислить 3 2015 с абсолютной погрешностью = 10−5. Это можно представить как нахождения нуля функции( ) = 2015 − 3. Применяя метод бисекции на отрезке [0, 100], получаем следующий результат:

144

a = 0.0000000, b = 100.0000000, c = 50.0000000, f(c) = -122985.0000000 a = 0.0000000, b = 50.0000000, c = 25.0000000, f(c) = -13610.0000000

a = 0.0000000, b = 25.0000000, c = 12.5000000, f(c) = 61.8750000

a = 12.5000000, b = 25.0000000, c = 18.7500000, f(c) = -4576.7968750 a = 12.5000000, b = 18.7500000, c = 15.6250000, f(c) = -1799.6972656 a = 12.5000000, b = 15.6250000, c = 14.0625000, f(c) = -765.9143066 a = 12.5000000, b = 14.0625000, c = 13.2812500, f(c) = -327.7009583 a = 12.5000000, b = 13.2812500, c = 12.8906250, f(c) = -127.0121193 a = 12.5000000, b = 12.8906250, c = 12.6953125, f(c) = -31.1156964 a = 12.5000000, b = 12.6953125, c = 12.5976562, f(c) = 15.7400736

a = 12.5976562, b = 12.6953125, c = 12.6464844, f(c) = -7.5973567 a = 12.5976562, b = 12.6464844, c = 12.6220703, f(c) = 4.0939285 a = 12.6220703, b = 12.6464844, c = 12.6342773, f(c) = -1.7460661 a = 12.6220703, b = 12.6342773, c = 12.6281738, f(c) = 1.1753425 a = 12.6281738, b = 12.6342773, c = 12.6312256, f(c) = -0.2850089 a = 12.6281738, b = 12.6312256, c = 12.6296997, f(c) = 0.4452550 a = 12.6296997, b = 12.6312256, c = 12.6304626, f(c) = 0.0801451 a = 12.6304626, b = 12.6312256, c = 12.6308441, f(c) = -0.1024264 a = 12.6304626, b = 12.6308441, c = 12.6306534, f(c) = -0.0111393 a = 12.6304626, b = 12.6306534, c = 12.6305580, f(c) = 0.0345033 a = 12.6305580, b = 12.6306534, c = 12.6306057, f(c) = 0.0116821 a = 12.6306057, b = 12.6306534, c = 12.6306295, f(c) = 0.0002714 a = 12.6306295, b = 12.6306534, c = 12.6306415, f(c) = -0.0054339 a = 12.6306295, b = 12.6306415, c = 12.6306355, f(c) = -0.0025813 a = 12.6306295, b = 12.6306355, c = 12.6306325, f(c) = -0.0011549 a = 12.6306295, b = 12.6306325, c = 12.6306310, f(c) = -0.0004417 a = 12.6306295, b = 12.6306310, c = 12.6306303, f(c) = -0.0000852 a = 12.6306295, b = 12.6306303, c = 12.6306299, f(c) = 0.0000931 a = 12.6306299, b = 12.6306303, c = 12.6306301, f(c) = 0.0000040

c = 12.6306301

Можно убедиться, что значение получено верно. Оценим скорость сходимости данного метода. На каждом шаге отрезок [ , ] уменьшается вдвое. Пусть — середина отрезка [ , ], полученная на -ом шаге, а — искомое решение. Тогда на -ом шаге, очевидно, выполняется:

	−		| 6	−	.
|				2

Из этого неравенства можно найти необходимое количество шагов для достижения тре-

буемой погрешности : = log2 ( − ). Таким образом, метод сходится относительно быстро

(хотя и не очень) и позволяет найти корень с любой наперед заданной точностью. Он нам пригодится в дальнейшем.

2 Нахождение всех корней полинома

Предположим, что имеется многочлен ( )	= 0 + 1 + 2 2 + · · · + над полем R,
и необходимо найти все такие вещественные ,	что ( ) = 0. По теореме Абеля-Руффини

данная задача неразрешима в радикалах при > 4. Однако искать корни многочленов приходится во многих областях, поэтому существует великое множество предназначенных для этого

145

численных методов. Мы изучим не самый эффективный, но наиболее простой метод решения этой задачи. Как и раньше, мы предполагаем, что задано некоторое и мы хотим найти такие

: | ( )| < .

Предположим, что есть некоторые 2 точки и , причем полином имеет разные знаки в этих точках. Тогда мы можем воспользоваться уже описанным методом бисекции, чтобы найти какой-либо корень. Однако нам необходимо найти все корни. Заметим, что если функция ( ) строго монотонна на отрезке [ , ], то на этом отрезке существует единственный корень. Действительно, если предположить, что существуют две такие точки 1 < 2 : ( 1) = ( 2) = 0, то в силу монотонности ( 1, 2) 0 < ( ) < 0, что невозможно. В данном случае мы можем с полной уверенностью требовать именно строгую монотонность, т.к. если предположить, что имеется отрезок [ ′, ′], такой, что полином равен некоторому во всех точках этого отрезка, то выходит, что у полинома ( ) − корней бесконечно много, что противоречит основной теореме алгебры.

Для того, чтобы на отрезке [ , ] функция была бы монотонна, очевидно, необходимо и достаточно того, чтобы и были экстремумами, причем @ : < < и — экстремум. Значит, мы свели исходную задачу к нахождению всех экстремумов ( ). По теореме Ферма, если 0 — экстремум функции ( ), то ′( 0) = 0, причем обратное, вообще говоря, неверно. В нашем случае этому свойству будут также удовлетворять все точки перегиба.

Предположим, что мы нашли все такие (в количестве штук), что ′ ( ) = 0. Расположим их так, что = 1, − 1 < +1. Заметим, что точки перегиба не меняют характера монотонности. Иными словами, если есть экстремумы и и точка перегиба между ними, то мы ничего не потеряем, если применим бисекцию на отрезке [ , ], а затем на [ , ], т.к. если функция монотонна на [ , ], то она монотонна и на [ , ], и на [ , ]. Таким образом, чтобы найти все корни, нужно просто применить метод бисекции на тех интервалах (−∞, 1], [ 1, 2], . . . , [ −1, ], [ , +∞), где полином имеет разные знаки в границах. Некоторую сложность может представлять нахождение «конкретных» значений ±∞, но обычно корни ищутся на некотором отрезке [ , ]. Если же это не так, то соответствующие границы определяются вручную (чаще всего полиномы растут достаточно быстро, поэтому с этим проблем нет). Существуют и конструктивные методы определения границ корней, которые читатель может найти самостоятельно.

Соответственно, осталось лишь научиться находить все такие . Как ни странно, это делается очень просто. Заметим, что ′ ( ) — многочлен уже ( − 1)-ой степени. Нам необходимо найти все его нули. Таким образом, получили простую рекурсию:

Пусть функция ( ) возвращает список всех нулей полинома .

База: = 1. В этом случае мы либо тривиально находим единственный корень и возвращаем его, либо возвращаем пустой список.

В противном случае продифференцируем ′ ( ) = −1( ). Пусть теперь список ̂ =

( −1) длины .

Применим метод бисекции с нашим последовательно ко всем интервалам (−∞, ̂1], [̂1, ̂2], . . . , [̂ , +∞). Пусть найденные решения 1, . . . , ′ (в порядке обработки интервалов). Вернем в качестве результата список ( 1, . . . , ′).

146

Например,

возьмем

многочлен

3( )

=

3

2

После взятия производ-

2 −

3 − 72 + 35.

2

− 6 − 72. Корнями

ной получаем 2( ) = 6

данного многочлена являются 1 = −3, 2 = 4.

Достаточно применить 3 раза метод бисекции,

чтобы найти все 3 корня.

Теперь оценим время работы данного мето-

да. Очевидно, что время работы каждого вызо-

ва метода бисекции можно оценить сверху как

(

(

−

) )

O

log

2

, т.к. для вычисления значения полинома в точке требуется O( ) операций

(по схеме Горнера).

Тогда если ( ) — время работы функции для полинома , то

−

·

(

2 (

) )

( ) = (

1) +

O

log

−

,

(

2

(

) )

( ) = O

3

log

−

.

3Метод Ньютона

На практике описанный выше метод бисекции применяется достаточно редко ввиду его малой скорости сходимости. Если необходимо вычислить значение какой-либо элементарной функции, вроде корня, многие математические пакеты, библиотеки языков программирования или калькуляторы обычно используют либо разложение в ряд Тейлора с последующим взятием первых нескольких членов, либо производят вычисление по методу Ньютона. Несмотря на свой почтенный возраст (что можно понять по названию), метод касательных (второе название) до сих пор используется там, где необходима высокая скорость вычислений (или их тривиальность). Яркий пример: алгоритм быстрого вычисления обратного квадратного корня

√

1 . Впрочем, применять метод Ньютона нужно с большой осторожностью, подробнее о чем

будет сказано ниже.

Пусть есть функция ( ): [ , ] →R, дифференцируемая на отрезке [ , ], и нужно найти произвольное решение уравнения ( ) = 0. Сначала дадим геометрическое описание метода. На рисунке справа приведен пример функции (выделена синим), нуль которой нужно найти. Пусть — некоторая точка, в которой значение функции отлично от нуля. Проведем касательную в этой точке и найдем ее пересечение с осью . Обозначим его как +1. Как видим, новая вычисленная точка лежит ближе к искомому корню. Как мы уже знаем,

′( ) = tg =				=		( )	=	− ( )	,
						− +1
								+1 −
+1 = −	( )				.
	′( )

Последняя формула и задает последовательность, по которой мы вычисляем приближения. Таким образом, метод Ньютона заключается в задании начального приближения 0 и последовательном вычислении до тех пор, пока не станет верным неравенство | ( )| < (или

147

| +1 − | < , что лучше, т.к. решения может и не быть). Например, найдем решение уравнения cos = 3 с абсолютной погрешностью = 10−12. Необходимо найти нуль функции( ) = cos − 3. Найдем производную ′( ) = − sin − 3 2. В качестве начального приближения возьмем 0 = 0.5, т.к. искомый корень, очевидно, не превосходит 1. Получим следующую последовательность приближений:

0 = 0.500000000000

1 = 1.112141637097

2 = 0.909672693737

3 = 0.867263818209

4 = 0.865477135298

5 = 0.865474033111

6 = 0.865474033102

Итоговый результат получен достаточно быстро. Число точных разрядов примерно удваивается с каждым шагом, поэтому метод обладает квадратичной сходимостью. На рисунке справа приведено построение нескольких первых приближений.

Теперь рассмотрим уже знакомый пример вычис-

√

ления 3 2015 с точностью 10−6. Снова представим это как нахождение нуля функции ( ) = 2015 − 3. В качестве начального приближения возьмем, например, точку 0 = 50. По формуле получаем

	+1 = +	2015 − 3	.
		3 2
	Тогда имеем следующие приближения:
0 = 50.0000000
1 = 33.6020000
2 = 22.9962054
3 = 16.6009139
4	= 13.5044682
5	= 12.6859542
6	= 12.6308710
7	= 12.6306301

Итоговое решение найдено почти в 3 раза быстрее, чем при использовании метода бисекции, при том, что поиск проводился при меньшей погрешности! Однако продемонстрируем теперь пример, в котором возможно столкнуться с проблемой при применении метода Ньютона. Рассмотрим функцию ( ) = 3 − 2 + 2. Теперь возьмем в качестве начального приближения0 = 0. Следующая вычисленная точка будет равна1 = 1, после чего 2 = 0. Таким образом метод никогда не сойдется, хотя решение и существует.

148