Схема Горнера и метод Ньютона

Если известны корни характеристического уравнения n порядка (1.10), то его можно представить в виде:

, (1.11)

где x₀ , x₁ , ¼, x_n – корни уравнения. Исходный полином (1.10) разложили на множители. Соотношение между коэффициентами и корнями уравнения можно получить перемножив левую часть уравнения (1.11). Этот алгоритм существует.

Аналитическое решение известно лишь для уравнений до 4-го порядка. Численные методы позволяют рассчитать корни характеристического уравнения высокого порядка. Знание корней необходимо для получения аналитического решения обыкновенных дифференциальных уравнений и для анализа устойчивости систем.

Критерии устойчивости разрабатывались по той причине, что они позволяли судить об устойчивости системы, не вычисляя корней. А вычислить корни без ЭВМ было весьма трудно.

Рассмотрим наиболее распространенные методы численного решения алгебраических уравнений. Для этого введем в рассмотрение схему Горнера, которая позволяет более рационально рассчитывать полиномы вида (1.10):

.

Запишем правую часть функции f(x) в виде:

(1.11)

Если раскрыть скобки и приравнять коэффициенты при одинаковых степенях x в формулах (1.10) и (1.11), то получим соотношение между коэффициентами a и b:

(1.12)

коэффициенты a известны, поэтому перепишем соотношения (1.12) в виде:

(1.13)

где j=1, 2, …, n. Если x₀ корень уравнения, то мы получили соотношение (1.13) для коэффициентов уравнения порядка. Используем его для записи уравнения 5-го порядка: (1.14)

Получен алгоритм расчета полинома по схеме Горнера - последняя строка (1.14). Если x₀ корень уравнения (1.10), то коэффициент b_n = 0.

Поделим обе части (1.11) на x-x₀:

Если x₀ корень уравнения (1.10), то деление полиномов происходит без остатка (теорема Безу). Остаток (для любого x₀) определяется по схеме Горнера (последняя строка (1.14)). Для уравнения n порядка вычисление по схеме Горнера выполняются по алгоритму:

(1.15)

Обозначим f(x)=y и запишем каждую пару слагаемых, начиная с последней, следующим образом:

присвоим y = a₀ и в цикле считаем

y=a(k)+y*x, где k=1, 2, 3, …, n. (1.16)

При вычислении по схеме Горнера требуется 2n арифметических действий, а по формуле (1.10) .

Реализуем алгоритм на МATLAB.

function [f]=gornR

a = [1 3 3 1]; n = 3; y = a(1); x = -1;

for k=2: n+1

y=a(k)+y*x;

end

Конечно, алгоритм справедлив лишь для вещественных корней. Его можно применять для алгебраического уравнения, если n нечетное. Если корень определен, то далее используется алгоритм понижения порядка (1.13):

function [a, b]=ponporR

n = 4; a = [1 4 6 4 1]; x = -1;

b(1) = a(1);

for k=1: n-1

b(k+1) = a(k+1) + x*b(k);

end

n=n-1;

for k=1: n+1

a(k)=b(k);

end

Т.е. новый вектор a получендля характеристического уравнения степени .

Определить вещественный корень можно методом Ньютона. Алгоритм метода легко получить из графика

рис. 1.1

Задаем начальное приближение x₀, определяем f(x₀) и производную . Пересечение касательной с осью абсцисс позволяет определить новое приближение x₁.

Запишем приближенное значение производной в точке x₀

, откуда

В общем виде алгоритм запишется:

i - текущая итерация.

Процесс итераций по x следует повторять до достижения величиной заданной точности ε. Предыдущее и последующее значения x будут совпадать с точностью до ε.

Кроме значения функции f(x_i) требуется знать производную в текущей точке.

Получим ее, записав уравнение (1.11) в виде: где

Продифференцируем f(x), получим:

Пусть х=х_о, тогда, т.е. полином, пониженного на единицу порядка, является производной от исходного полинома.

Если использовать рекуррентное соотношение (1.16), то, продифференцировав его, получим следующий алгоритм:

присвоим

и в цикле считаем (1.18)

Объединив алгоритмы (1.16)-(1.18), реализуем их на МATLAB:

function [x]=newtonR;

global f x;

xа =-1.2; eps=0.0001;

for m=1: 20

[f]=gornerR;

x=x-f;

if (abs(f)< eps), break, end

end

function [f]=gornerR;

global f x

n=3; a=[1 3 3 1]; y=a(1); y1=0;

for k=2: n+1

y1=y+x*y1;

y=a(k)+x*y;

end

f=y./y1;

При обращении из МATLAB [x]=newtonR получим значение корня x = -1.0001.

В m-файле newtonR имеется две возможности остановить итерации: задавая параметр цикла 20, либо по условному оператору if. Количество итераций можно определить более простыми методами. Попробуйте это сделать самостоятельно.

Если все корни вещественные, то, использовав процедуру понижения порядка, можно определить новый корень и т.д. до тех пор пока порядок уравнения не будет равным 2. Решение квадратного уравнения позволит определить два оставшихся корня.

Весьма полезным для работы с полиномами является следующий m. файл:

function [s]=SumPoly(p, q)

%Определяется сумма векторов p и q. Они должны быть заданы

%Определяется длина векторов р и q и выбирается большая из них

max1=max(length(p), length(q));

%Фиксирутся длина нулевых векторов p1 и q1

p1=zeros(1, max1);

q1=zeros(1, max1);

%Формируются векторы p1 и q1 равной длины

p1(max1-length(p)+1: max1)=p;

q1(max1-length(q)+1: max1)=q;

s=p1+q1;

Задаем исходные полиномы произвольного порядка p и q: p=[1 1 1 1 1]; q=[2 2 2 2 2 2];

s – алгебраическая суммаполиномов: s =[2 3 3 3 3 3]. При составлении собственных программ преобразования передаточных функций приемы программирования, использованные в данном m.файле, могут оказаться очень полезными.

Умножение полиномов происходит с помощью оператора [a]= conv(p, q). Результат a – вектор коэффициентов полинома произведения.

Деление полиномов с помощью оператора [a, b]=deconv(p, q). Целая часть – а, коэффициенты полинома числителя остатка деления – b.

Коэффициенты полиномов можно определить по корням уравнения с помощью оператора poly(x).

Для счета функций можно использовать стандартные операторы eval и feval. Например, eval(‘x.^3+x.^2+x+1’) и feval(‘func’, 1). В первом случае предварительно следует задать x = 1, а во втором запасти в func.m файле необходимую функцию и задать значение аргумента. Функция может зависеть от многих переменных.

Задание для самостоятельной работы:

Реализовать на МATLAB решение алгебраических уравнений

а) ,

б) ,

в) , т.е. определить все корни.

г) Дать пояснения для операторов SumPoly(p, q).m файла.

д) Поработать с операторами conv, deconv, poly, eval, feval.

Рассмотрим общий случай, когда корни уравнения (1.10) комплексные.

Обозначим и рассмотрим применение схемы Горнера к расчету полинома f(x). Подставив комплексный корень в f(x), мы получим комплексную функцию f(x₁, x₂). Вещественную часть функции обозначим y₁, а мнимую y₂, т.е. y=f(x)=y₁+jy₂ .

Запишем схему Горнера, используя (1.16) следующим образом:

(1.19)

Рассмотрим алгоритм расчета производной для комплексного x с использованием схемы Горнера. Обозначим вещественное значение производной F₁, а мнимое значение F₂. Дифференцируя соотношения (1.19), получим рекуррентные соотношения:

(1.20)

Алгоритмы расчета y=f(x) и для комплексных корней получены. Используем их для алгоритма Ньютона (1.17). Если отношение к привести к алгебраической форме записи, то получим два уравнения:

(1.21)

которые позволяют определить вещественную и мнимую части искомого корня . Итерации обозначены индексами i.

В представленных далее m.файлах использованы полученные алгоритмы, но с помощью их решение можно получить лишь для одного корня. Необходимо самостоятельно реализовать алгоритмы для определения всех корней алгебраического уравнения для приведенного далее примера. Требуется также объяснить причины возникающих при этом трудностей.

function [y, f]=gornerC;

%Расчет полинома по схеме Горнера для комплексных корней

% %y(1)-вещественная часть полинома, y(2)-мнимая часть полинома,

%f(1)-вещественная часть производной от полинома,

%f(2)-мнимая часть производной от полинома.

global F x;

n=4; a=[1 5 10 10 4];

y(1)=a(1); y(2)=0; f(1)=0; f(2)=0;

for k=2: n+1

R(1)=y(1)+x(1)*f(1)-x(2)*f(2);

f(2)=y(2)+x(1)*f(2)+x(2)*f(1); f(1)=R(1);

R1(1)=a(k)+x(1)*y(1)-x(2)*y(2);

y(2)=x(1)*y(2)+x(2)*y(1); y(1)=R1(1);

end;

%Вещественная часть отношения функции к производной от нее в алгоритме

%Ньютона F(1) и F(2)-мнимая часть.

F(1)=y(1)*f(1)+y(2)*f(2);

F(2)=y(2)*f(1)-y(1)*f(2);

%Знаменатель

Zn(1)=f(1)^2+f(2)^2;

%Отношение

F(1)=F(1)/Zn(1); F(2)=F(2)/Zn(1);

function [x]=newtonC;

%Расчет комплексного сопряженного корня методом Ньютона

%Обращение из командной строки [x]=newtonC

global F x;

x(1)=-1.2; x(2)=1.0; eps=0.001;

; %Для цикла с предусловием необходимо задать заведомо неприемлемые значения

F(1)=10; F(2)=10;

%Далее следуют два варианта заданной точности счета

while sqrt(F(1)*F(1)+F(2)*F(2))> eps % abs(F(1)+F(2))> eps

[F]=gornerC;

x(1)=x(1)-F(1);

x(2)=x(2)-F(2);

end

Идентификаторы F1 и F2 использовали в файле gornerC. m дважды.

Значения этих параметров из файла gornerC.m передаются в файл newtonC.m (1.21).

Для комплексных сопряженных корней х=х₁± jх₂ запишем квадратное уравнение:

x² + рx + q = 0,

где p = - 2x₁, q =x₁² +x₂².

Поэтому исходное уравнение (1.10) можно записать в виде:

y(x) = (x² + px + q)(b_n-2 + b_n-3x+…+b₁x^n-3 + b₀x^n-2), (1.22)

если x=x₁ ± jx₂ – корни исходного полинома y(x).

Раскрывая скобки и приравнивая коэффициенты при одинаковых степенях x у исходного (1.10) и правого полиномов (1.22) получим:

b₀ = a₀,

b₁ = a₁ - pb₀

b₂ = a₂ - pb₁ - qb₀

……………….

b_n-2 = a_n-2 - pb_n-3-qb_n-4

Коэффициенты определяются рекуррентно. Запишем полученный алгоритм в виде процедуры:

function [a, b]=ponporC;

%Заданы размер вектора коэффициентов полинома n и значение комплексного корня

% x(1)+-j x(2)

%Получили вектор коэффициентов b для полинома степени n-2 и запасли его в вектор a

%%Обращение из командной строки [a, b]=ponporC

n=5; a=[1 4 7 6 2]; b=zeros(1, n);

x(1)=-1; x(2)=1;

b(1)=a(1); p=-2*x(1); q=x(1)^2+x(2)^2;

b(2)=a(2)-p*b(1);

for k=3: n-2

b(k)=a(k)-p*b(k-1)-q*b(k-2);

end;

for k=1: n-2

a(k)=b(k);

end

for k=n: -1: n-1

a(k)=[];

end

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

Сходимость метода Ньютона зависит от выбора начального приближения для искомого корня. Точность определения корней ε задается. Обычно используется неравенство . Точность может быть задана и в относительных величинах. Можно задавать точность определения искомого корня

Недостатки метода Ньютона.

Наличие в знаменателе алгоритма Ньютона производной приводит к тому, что для кратных корней одновременно и

y = f(x), и равны 0. Аналогичная ситуация может возникнуть для пары близко расположенных корней.

Возможно также зацикливание итераций. Эти варианты показаны на рис 1.2 а) и б).

а) б)

рис.1.2

Задавая начальное приближение по x в точке a и, в соответствии с алгоритмом Ньютона, двигаясь по касательной к кривой, попадаем в точку b. Откуда, по касательной к кривой, вновь приходим в точку a. Происходит «зацикливание» алгоритма расчета.

Поэтому алгоритм Ньютона работает не всегда. В таких случаях говорят, что метод не сходится.

Задание для самостоятельной работы:

Реализовать на МATLAB решение алгебраических уравнений

а) ,

б) ,

в) ,

г) ,

д) , т.е. определить все корни.

Рассмотрим метод решения алгебраических уравнений, который хорошо работает при наличии кратных корней.