МИНИСТЕРСТВООБРАЗОВАНИЯ И НАУКИ УКРАИНЫ
СУМСКИЙГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ
КАФЕДРАИНФОРМАТИКИ
Контрольнаяработа
ПОЧИСЛЕННЫМ МЕТОДАМ
натему:
“Методы нахождения корнейполиномов”
Сумы- 2007 г.
/>/>/>/>План1 Нахождение корней уравнений (EquationSection 1)
2 Схема Горнера3Функции произвольного вида
4 Нахождение корней полиномов
Список используемыхинформационных источников
1 Нахождениекорней уравнений (Equation Section 1)
Одним из наиболее распространенных методов поиска корней уравненийявляется метод Ньютона и его модификации. Пусть требуется решить уравнение/>. Будем считать,что xявляется решением уравнения. Разложим функцию f(x) в ряд в точке x0 близкой к точке x и ограничимся только первымидвумя членами разложения.
/>
Поскольку x – корень уравнения, то />. Следовательно,
/>
Таким образом, если нам известно приближенное значение корня уравнения,то полученное уравнение позволяет его уточнить. Понятно, что процесс уточненияможно повторять многократно, до тех пор, пока значение функции не будутотличаться от нуля на величину меньшую, чем заданная точность поиска. Очередноеk-е приближение находитсяпо формуле
/>
Ограничившись в разложении только первыми двумя членами, мы фактическизаменили функцию f(x)на прямую линию, касательную в точке x0, поэтому метод Ньютона еще называют методомкасательных. Далеко не всегда бывает удобно находить аналитическое выражениедля производной функции. Однако, в этом и нет особой необходимости: посколькуна каждом шаге мы получаем приближенное значение корня, можно для его вычисленияиспользовать приближенное значение производной.
/>
В качестве малой величины />можно взять, например, заданную точностьвычислений />,тогда расчетная формула примет вид
/>(1.1)
Сдругой стороны, для вычисления производной можно воспользоваться значениямифункции, полученными на двух предыдущих шагах,
/>(1.2)
Втаком виде метод называется методом секущих (secant method). При этом, однако,возникает проблема с вычислением первого приближения. Обычно полагают, что />, то естьпервый шаг вычислений проводится с использованием формулы (1.1), а всепоследующие – с использованием формулы (1.2). Именно эта вычислительнаясхема реализована в пакете Mathcad. Используя метод секущих, мы не можемгарантировать, что корень находится между двумя последними приближениями.Можно, однако, для вычисления очередного приближения использовать границы интервала,на котором функция меняет знак. Такой метод называется методом хорд (false position method).
Идеяметода секущих развивается в методе Мюллера. Однако в этом методе длянахождения очередного приближения используются три предыдущие точки. Инымисловами, метод использует не линейную, а квадратичную интерполяцию функции.Расчетные формулы метода следующие[1]:
/>(1.3)
/>(1.4)
Знакперед корнем выбирается так, чтобы абсолютное значение знаменателя быломаксимальным.
Посколькупоиск корня заканчивается, когда выполнится условие/>, то возможно появление ложныхкорней. Например, для уравнения /> ложный корень />появится в том случае,если точность поиска задана меньше, чем 0,0001. Увеличивая точность поиска,можно избавиться от ложных корней. Однако не для всех уравнений такой подходработает. Например, для уравнения />, которое, очевидно, не имеетдействительных корней, для любой, сколь угодно малой точности найдется значениеx, удовлетворяющеекритерию окончания поиска. Приведенные примеры показывают, что к результатамкомпьютерных вычислений всегда нужно относиться критически, анализировать их направдоподобность. Чтобы избежать «подводных камней» при использованиилюбого стандартного пакета, реализующего численные методы, нужно иметь хотя быминимальное представление о том, какой именно численный метод реализован длярешения той или иной задачи.
В том случае, когда известен интервал, на котором расположенкорень, можно воспользоваться иными методами нахождения решения уравнения.
В методе Риддера (Ridder’s method) вычисляют значение функции в середине интервала/>. Затемищут экспоненциальную функцию такую, что /> Затем применяют метод хорд,используя значения />. Очередное значение вычисляют поформуле
/>(1.5)
Метод Брента (Brent method) соединяет быстроту метода Риддера игарантированную сходимость метода деления отрезка пополам. Метод использует обратнуюквадратичную интерполяцию, то есть ищет x как квадратичную функциюy. На каждом шагепроверяется локализация корня. Формулы метода достаточно громоздки и мы небудем их приводить.
Особые методы применяют для поиска корней полинома. В этом случаемогут быть найдены все корни. После того как один из корней полинома найден,степень полинома может быть понижена, после чего поиск корня повторяется.
Метод Лобачевского, метод приближённого (численного) решенияалгебраических уравнений, найденный независимо друг от друга бельгийскимматематиком Ж. Данделеном, русским математиком Н. И. Лобачевским (в 1834 внаиболее совершенной форме) и швейцарским математиком К. Греффе. Суть Л. м.состоит в построении уравнения f1(x) = 0, корни которого являются квадратами корней исходногоуравнения f(x) = 0. Затем строятуравнение f2(x) = 0, корнями которогоявляются квадраты корней уравнения f1(x) = 0. Повторяя этот процесс несколько раз, получаютуравнение, корни которого сильно разделены. В случае если все корни исходного уравнениядействительны и различны по абсолютной величине, имеются простые вычислительныесхемы Л. м. для нахождения приближённых значений корней. В случае равных поабсолютной величине корней, а также комплексных корней вычислительные схемы Л.м. очень сложны.
Метод Лагерра (Laguerre’s method) основывается на следующих соотношениях дляполиномов
/>
Предполагают, что корень x1 находится на расстоянии a от текущего приближения,в то время как все другие корни находятся на расстоянии b: />. Тогда
/>,
откуда
/>
Знак перед корнем выбирают с таким расчетом, чтобы получить наибольшеезначение знаменателя.
Еще один метод, который применяют для поиска корней полиномов, –метод сопровождающей матрицы (companion matrix). Можно доказать, чтоматрица
/>,
называемая сопровождающейматрицей для полинома />, имеет собственные значенияравные корням полинома. Напомним, что собственными значениями матрицыназываются такие числа , для которых выполняется равенство /> или />. Существуют весьмаэффективные методы поиска собственных значений, о некоторых из них мы будемговорить далее. Таким образом, задачу поиска корней полинома можно свести кзадаче поиска собственных значений сопровождающей матрицы.
2Схема Горнера
Вычисление по схеме Горнера оказывается более эффективным, причемоно не очень усложняется. Эта схема основывается на следующем представлениимногочлена:
p(x)= ((… ((anx + an-1)x + an-2)x +… + a2)x + a1)x + a0.
Займемся общим многочленом вида:
p(x) = anxn + an-1xn-1 + an-2xn-2 +… + a2x2 + a1x + a0.
Мы будем предполагать, что все коэффициенты an, ..., a0 известны,постоянны и записаны в массив. Это означает, что единственным входным даннымдля вычисления многочлена служит значение x, а результатом программы должнобыть значение многочлена в точке x.
Стандартный алгоритм вычисления прямолинеен:
Evaluate(x)
xточка, в которой вычисляется значение многочлена
result = a[0] + a[1]*x
xPower = x
for i = 2 to n do
xPower = xPower*x
result = result + a[i]*xPower
end for
return result
Этот алгоритм совершенно прозрачен и его анализ очевиден. В циклеfor содержится два умножения, которые выполняются n — 1 раз. Кроме того, одноумножение выполняется перед циклом, поэтому общее число умножений равно 2n — 1.В цикле выполняется также одно сложение, и одно сложение выполняется передциклом, поэтому общее число сложений равно n.
Вы можете легко проверить, что это представление задает тот жемногочлен, что и выше. Соответствующий алгоритм выглядит так:
HornersMethod(x)
xточка, в которой вычисляется значение многочлена
for i = n — 1 down to 0 do
result = result*x
result = result + a[i]
end for
return result
Цикл выполняется n раз, причем внутри цикла есть одно умножение иодно сложение. Поэтому вычисление по схеме Горнера требует n умножениё и nсложений — двукратное уменьшение числа умножений по сравнению со стандартнымалгоритмом.
Предварительная обработка коэффициентов
Результат можно даже улучшить, если обработать коэффициентымногочлена до начала работы алгоритма. Основная идея состоит в том, чтобывыразить многочлен через многочлены меньшей степени. Например, для вычислениязначения x256 можно воспользоваться таким же циклом, как и в функции Evaluate вначале этой статьи; в результате будет выполнено 255 умножений. Альтернативныйподход состоит в том, чтобы положить result= x*x, а затем семь раз выполнить операцию result = result*result. Послепервого выполнения переменная result будет содержать x4, после второго x8,после третьего x16, после четвертого x32, после пятого x64, после шестого x128,и после седьмого x256.
Для того, чтобы метод предварительной обработки коэффициентовработал, нужно, чтобы многочлен был унимодальным (то есть старший коэффициентan должен равняться 1), а степень многочлена должна быть на единицу меньшенекоторой степени двойки (n = 2k — 1 для некоторого k > 1). В таком случаемногочлен можно представить в виде:
p(x)= (xj + b)q(x) + r(x), где j = 2k-1. (1)
В обоих многочленах q и r будет вдвое меньше членов, чем в p. Дляполучения нужного результата мы вычислим по отдельности q(x) и r(x), а затемсделаем одно дополнительное умножение и два сложения. Если при этом правильновыбрать значение b, то оба многочлена q и r оказываются унимодальными, причемстепень каждого из них позволяет применить к каждому из них ту же самуюпроцедуру. Мы увидим, что ее последовательное применение позволяет сэкономитьвычисления.
Вместо того, чтобы вести речь о многочленах общего вида, обратимсяк примеру. Рассмотрим многочлен:
p(x)= x7 + 4x6 — 8x4 + 6x3 + 9x2 + 2x — 3.
Определим сначала множитель xj + b в уравнении (1). Степеньмногочлена p равна 7, то есть 23 — 1, поэтому k = 3. Отсюда вытекает, что j =22 = 4. Выберем значение b таким образом, чтобы оба многочлена q и r былиунимодальными. Для этого нужно посмотреть на коэффициенты aj-1 в p и положить b= aj-1 — 1. В нашем случае это означает, что b = a3 — 1 = 5. Теперь мы хотимнайти значения q и r, удовлетворяющие уравнению:
x7 + 4x6 — 8x4 + 6x3 + 9x2 + 2x — 3 = (x4 + 5)q(x) + r(x).
Многочлены q и r совпадают соответственно с частным и остатком отделения p на x4 + 5. Деление с остатком дает:
p(x)= (x4 + 5)(x3 + 4x2 + 0x + 8) + (x3 — 11x2 + 2x — 37).
На следующем шаге мы можем применить ту же самую процедуру ккаждому из многочленов q и r:
q(x) = (x2 — 1)(x + 4) + (x + 12), r(x) = (x2 + 1)(x — 11) + (x — 26).
В результате получаем:
p(x)= (x4 + 5)((x2 — 1)(x + 4) + (x + 12)) + ((x2 + 1)(x — 11) + (x — 26)).
Посмотрев на этот многочлен, мы видим, что для вычисления x2требуется одно умножение; еще одно умножение необходимо для подсчета x4 =x2*x2. Помимо этих двух умножений в вычислении правой части равенства участвуютеще три умножения. Кроме того, выполняется 10 операций сложения. В таблице 1приведены сравнительные результаты анализа этого метода и других методоввычисления. Экономия не выглядит значительной, однако это объясняется тем, чтомы рассматриваем лишь частный случай. Общую формулу для сложности можновывести, внимательно изучив процедуру. Заметим прежде всего, что в равенстве(1) участвуют одно умножение и два сложения. Поэтому для числа умножений M =M(k) и числа сложений A = A(k) мы получаем следующий набор рекуррентныхсоотношений: M(1) = 0 A(1) = 0 M(k) = 2M(k — 1) + 1 при k > 1 A(k) = 2A(k — 1) + 2 при k > 1.
Таблица1. Число операций при вычислении значения многочлена степени 7
Способ
Умножения
Сложения Стандартный 13 7 Схема Горнера 7 7 Предварительная обработка 5 10
Решая эти соотношения, мы заключаем, что число умноженийприблизительно равно N/2, а число сложений приблизительно равно (3N — 1)/2.Неучтенными, однако, остались умножения, необходимые для подсчетапоследовательности значений x2, x4, x8, ..., x2k-1; на это требуетсядополнительно k — 1 умножение. Поэтому общее число умножений приблизительноравно N/2 + log2N.
Таблица 2. Число операций при вычислениизначения многочлена степени N
Способ
Умножения
Сложения Стандартный 2N — 1 N Схема Горнера N N Предварительная обработка N/2 + log2N (3N — 1)/2
В таблице 2 приведены результаты сравнительного анализастандартного алгоритма, схемы Горнера и алгоритма с предварительной обработкойкоэффициентов. При сравнении последних двух алгоритмов видно, что нам удалосьсэкономить N/2 — log2N умножений, но за счет дополнительных (N — 1)/2 сложений.Во всех существующих вычислительных системах обмен умножений на сложения считаетсявыгодным, поэтому предварительная обработка коэффициентов повышаетэффективность./>/>3 Функции произвольного вида
Найдем нули функции /> на интервале x=[–2,7], используяMathcad
Изобразим сначала функцию на графике.
/>
/>
На заданном интервале функция три разаобращается в ноль. Определим нули функции, используя встроенную функциюroot(f(x),x). Первый аргумент – функция, нуль которой необходимо найти,второй – переменная, которую необходимо варьировать. (Вообще говоря,функция f может быть функцией многих переменных и необходимо указывать, покакой именно переменной мы ищем нуль функции.) Кроме того, необходимо задать начальноеприближение поиска. Точность вычислений задается встроенной переменной TOL. Поумолчанию ее значение равно 0,001. Это значение можно изменить либо через менюMath/Built–In Variables или непосредственно в тексте документа: />
Задаем начальное приближение: />
И вычисляем корень: />
Если требуется найти несколько корней, как внашей задаче, то имеет смысл определить новую функцию: />
Функция r(x) возвращает значение корняближайшее к x[2],то есть начальное приближение мы задаем через аргумент функции. Задаем векторначальных приближений x и находим соответствующие им корни X:
/> /> /> />
Для данного примера корни легко могут бытьнайдены аналитически. Они равны на заданном интервале /2/2и Полученный численныйрезультат с заданной точностью совпадает с точным решением.
Определение новой функции целесообразно и втом случае, когда мы хотим исследовать зависимость решения от параметра. Пустьфункция зависит от параметра a
/>
Первый аргумент функции z задает значениепараметра, второй – начальное приближение. Найдем корни уравнения призначениях параметра 1 и 2.
/> />
Если мы хотим получить комплексный корень,то начальное приближение следует задавать комплексным:
/>/>/>4 Нахождение корней полиномов
Для нахождения корней полиномов имеетсявстроенная функция polyroots(a). Аргументом функции является векторкоэффициентов полинома />, то есть для уравнения />вектор а имеетвид
/> /> />
Если в полиноме отсутствуют некоторыестепени, то на соответствующих местах следует писать 0. Пусть требуется найтикорни полинома />
/> /> />
Коэффициенты полинома могут быть икомплексными.
Списокиспользуемых информационных источников
1. Coppersmith D., Odlyzko A. M., Schroeppel R. Descrete logarithmsin GF(p) // Algorithmica. V. 1,1986. P. 1-15.
2. Lenstra A. K, Lenstra H. W. (jr.) The Development of the NumberField Siev. Lect. Notes in Math. V. 1554. Springer, 1993.
3. McCarthy D. P. “The optimal algorithm to evaluate xn usingelementary multiplication methods”, Math. Comp., vol. 31, no 137, 1977, pp. 251– 256.