ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2008, том 48, < 7, с. 1181-1186
УДК 519.615.5
ОПРЕДЕЛЕНИЕ КРАТНОСТИ КОРНЯ НЕЛИНЕЙНОГО АЛГЕБРАИЧЕСКОГО УРАВНЕНИЯ1)
© 2008 г. Н. Н. Калиткин*, И. П. Пошивайло**
(*125047 Москва, Миусская пл., 4а, ИММ РАН;
**124498 Зеленоград, пр-д 48065, Моск. гос. ин-т электронной техн.) e-mail: kalitkin@imamod.ru Поступила в редакцию 30.03.2007 г.
Переработанный вариант 29.01.2008 г.
Для нахождения корней нелинейного алгебраического уравнения наиболее часто используют метод Ньютона. Для расширения области сходимости метода Ньютона применяют одно обобщение, нередко называемое непрерывным аналогом метода Ньютона. Для классического и обобщенного методов Ньютона предложен эффективный метод нахождения корней с одновременным вычислением их кратности. При этом корни даже высокой кратности (до порядка 10) вычисляются с малой погрешностью. Метод проиллюстрирован численными примерами. Библ. 8. Табл. 2.
Ключевые слова: нелинейное алгебраическое уравнение, численное нахождение корней, определение кратности корней, обобщенный метод Ньютона.
1. ПРОБЛЕМА
Пусть /(х) - функция одного переменного, непрерывная вместе с некоторым числом своих производных. Для нахождения корней х* уравнения
/ (х) = 0 (1)
существует много неплохих методов. Они описаны в большом количестве учебников и монографий, например [1]-[3]. Однако до сих пор определение кратности найденного корня может вызвать серьезные затруднения. Рассмотрим их причины.
Очень надежный метод деления пополам в принципе не позволяет найти корни четной кратности. Корень нечетной кратности он находит, но без указания кратности. Если в этом случае пытаться исключать корни делением /(х) на х - х*, то после первого исключения останется корень четной кратности. Далее метод деления пополам перестает замечать этот корень и не позволяет найти его вторично и исключить.
Метод Ньютона в достаточно малой окрестности корня любой кратности сходится к нему. Поэтому после первого исключения найденного кратного корня он может снова к нему сходиться. Однако на практике мы находим не точное значение корня х*, а приближенное значение х* .
Поэтому второй раз мы находим корень с меньшей точностью. Для корней кратности порядка 5 и выше даже при расчетах на 64-разрядном компьютере может возникнуть значительная погрешность.
В прочих методах тоже возникают те или иные трудности. В литературе и программах для произвольной функции /(х) нам удалось найти только один нестрогий способ (см. [4]), позволяющий оценить кратность корня. О нем будет сказано далее.
Существует литература, специально посвященная вычислению корней многочлена (например, [5]). Это связано в основном с проблемой нахождения собственных значений матриц высокого порядка, которая сводится к вычислению корней характеристического полинома соответствующей степени. Для этой задачи разработаны некоторые специфические приемы. Однако
1 -'Работа выполнена при финансовой поддержке РФФИ (коды проектов 05-01-00152, 05-01-08006) и НШ-5772.2006.1.
они не переносятся на функции произвольного вида. На особенностях этой проблемы мы остановимся позднее.
Ниже предложено несложное дополнение к методу Ньютона и обобщенному методу Ньютона, позволяющее для функции произвольного вида находить корни высокой кратности с малой погрешностью и одновременно определять их кратность. Это дополнение пригодно для написания хороших стандартных программ.
2. МЕТОД НЬЮТОНА Классический метод Ньютона основан на итерационном процессе
x. + 1 = x. - f (x.) / f' (x.). (2)
Как известно, он сходится, если существует непрерывная функция f"(x), а начальное приближение выбрано в достаточно малой окрестности корня. Пусть существует непрерывная функция f(p + Х)(х), а x* естьр-кратный корень уравнения (1). Тогда в малой окрестности корня имеет место приближенная формула
f(x)~ a Ар + ЬАр + 1, А = x - x*. (3)
Подставляя (3) в (2) и вычитая из обеих частей равенства значение x*, получаем
4 + 1 = ^ As + O (А2).
Отсюда видно, что метод Ньютона для простого корня (р = 1) сходится квадратично, а для кратного корня (р > 2) - линейно. По скорости сходимости можно определить кратность корня.
Очевидно, что (р - 1)/р = A. + 1/Ai + O(As). Величины As, As + 1 нам не известны, поскольку в них входит неизвестное значение x*. Однако нетрудно показать, что таким же (с точностью до малых величин) будет отношение (xs + 1 - xs)/(xs - xs _ 1). Удобно ввести знаменатель линейной сходимости по формуле
qs = x±+:±Zx. « IL-l (4)
xs - xs -1 р
и величину
р- = 14; <5)
При сделанных предположениях рs —► р при s —►
Поэтому надо ввести в ньютоновские итерации расчет р. по (4) и (5). Если итерации сходятся, а р. стремится к целому числу р, то мы находим значение корня x* и его кратность. Корень определяется при этом с высокой точностью даже при большой кратности.
Заметим, что предложенный метод применим даже к функциям, имеющим вблизи корня вид f(x) ~ (x - x*)^ где р - не целое число. Сходимость значений р. к нецелому числу служит указанием на такой корень дробной кратности.
В [4] предложен следующий способ. Заменим (2) на следующий итерационный процесс:
x. + 1 = x. - mf (xs)/f'(xs), (6)
где m - натуральное число. Проведем расчеты с разными m, выбирая одно и то же нулевое приближение. То значение m, при котором число итераций окажется наименьшим, должно быть кратностью корня.
По существу, этот способ основан на тех же соображениях. Однако вместо анализа сходимости вблизи корня здесь берется число итераций, которое при далеком нулевом приближении может служить не слишком надежным критерием. Поэтому такой способ является гораздо менее надежным, чем предложенный здесь. Кроме того, этот способ не позволяет найти корни дробной кратности.
Таблица 1
е p 100 10-1/2 10-1 10-3/2 10-2 10-5/2 0
1 0 0 0 0 0 0 0
2 0.5000 0.5201 0.5265 0.5285 0.5291 0.5293 0.5294
3 0.6667 0.6851 0.6909 0.6927 0.6933 0.6935 0.6936
4 0.7500 0.7656 0.7705 0.7720 0.7725 0.7727 0.7728
5 0.8000 0.8133 0.8175 0.8188 0.8192 0.8193 0.8194
6 0.8333 0.8448 0.8485 0.8496 0.8500 0.8501 0.8501
7 0.8571 0.8673 0.8705 0.8715 0.8718 0.8719 0.8719
8 0.8750 0.8840 0.8869 0.8878 0.8881 0.8882 0.8882
9 0.8889 0.8970 0.8996 0.9004 0.9007 0.9008 0.9008
10 0.9000 0.9074 0.9098 0.9105 0.9107 0.9108 0.9108
11 0.9091 0.9159 0.9180 0.9187 0.9189 0.9190 0.9190
12 0.9167 0.9229 0.9249 0.9256 0.9258 0.9258 0.9259
13 0.9231 0.9289 0.9308 0.9313 0.9315 0.9316 0.9316
14 0.9286 0.9340 0.9357 0.9363 0.9365 0.9365 0.9365
15 0.9333 0.9384 0.9401 0.9406 0.9407 0.9408 0.9408
3. ОБОБЩЕННЫЙ МЕТОД НЬЮТОНА
Этот метод в непрерывной формулировке предложен в [6], а его численная дискретная реализация - в [7] и последующих работах тех же авторов. Это обобщение сохраняет квадратичную сходимость вблизи простого корня и расширяет область допустимых нулевых приближений. Вместо (2) записывается обобщенный итерационный процесс (см. [7])
Xs +1 = Xs - Ts9(xs), Ф( x) = f(x)/f (x), 0 <Ts < 1. (7)
Значения Ts выбирают так, чтобы вдали от корня они были небольшими, а при попадании в малую окрестность корня стремились бы к 1. Мы будем использовать так называемый оптимальный шаг:
2 2
Ts = f2( xs) + 9 f2 (xs - Ф ( xs) ), 0 <е< 1. (8)
f (xs) + f (xs - Ф(xs))
Здесь 9 - управляющий параметр метода. С ростом 9 величина Ts монотонно возрастает от значения
2
f (xs)
f2(xs) + f 2(xs - Ф(xs))
до т. = 1. Очевидно, что при б = 1 (т. = 1) формулы (7), (8) переходят в классический метод Ньютона (2).
Если т. < 1, то формулу (7) можно рассматривать как численную реализацию дифференциального уравнения из [6] простейшей схемой Эйлера. Это устанавливает некоторую связь между обобщенным методом Ньютона и методом из [6]. Реально такая ситуация может встретиться при малых значениях б вдали от корня. Однако при не очень малых т. такого соответствия уже не будет. Поэтому, хотя дифференциальный метод из [6] всегда сходится к корню при любом нулевом приближении, обобщенный метод Ньютона таким свойством не обладает. Можно рассчитывать лишь на расширение области сходимости по сравнению с классическим методом Ньютона.
В методе (7), (8) можно также определять кратность корня по скорости сходимости итераций. Нетрудно показать, что в малой окрестности корня величина (4) связана с кратностью корня следующим соотношением:
. 11 + 8 (1 - 1 / р)". (9) Р 1+ (1-1/р )2р
Явную зависимость р(д) из (9) получить нельзя. Поэтому целесообразен следующий способ. Табулируем зависимость д(р, 8) для серии значений 8 (табл. 1). Сопоставляя расчетные значения д. с числами в табл. 1 при используемом в расчете параметре 8, определяем кратность корня.
Таблица 2
Метод Ньютона
P классический обобщенный, б = 10 1/2
погрешность корня x* погрешность значения fx*) Peff погрешность корня x* погрешность значения fx*) peff
1 1 х 10-17 0 1.00 10-17 0 1
2 8 х 10-9 9.2 х 10-17 2.00 7 х 10-7 2.2 х 10-13 2
3 1 х 10-5 4.9 х 10-16 3.01 1 х 10-5 2.6 х 10-16 3
4 6 х 10-4 7.3 х 10-15 4.03 1 х 10-3 5.5 х 10-14 4
5 8 х 10-3 2.4 х 10-13 5.06 5 х 10-3 2.9 х 10-14 5
6 2 х 10-2 1.6 х 10-13 6.08 3 х 10-2 1.1 х 10-14 6
7 2 х 10-2 1.1 х 10-15 7.04 6 х 10-2 5.8 х 10-13 7
8 1 х 10-1 1.3 х 10-13 7.95 6 х 10-2 5.2 х 10-15 8
4. СТРАТЕГИЯ РАСЧЕТА
Практика показала, что целесообразно начинать расчет при параметре б = 1, т.е. классическим методом Ньютона. Если S итераций не обеспечивают сходимости, то параметр б рекомендуется уменьшать в 101/2 раз, и снова выполнять не более S итераций. При этом целесообразно выполнять второй итерационный процесс не с последнего найденного приближения, а с исходного нулевого: если итерации плохо сходятся, то последнему приближению рискованно доверять. Если следующих S итераций снова окажется недостаточно, процедуру уменьшения б в 101/2 раз повторяют. Заметим, что табл. 1 построена именно для такого уменьшения б. Если сходимости не будет даже при б = 0.001, задачу следует считать непосильной для программы.
Описанный способ тестировался на разнообразных функциях и сравнивался с рядом известных стандартных программ FSOLVE и FZERO из пак
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.