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

2.1 Предварительный анализ функции

Для начала рассмотрим все случаи поведения функции полинома 3-ей степени.

Р3(x)=a3x3 ∙ a2x2 ∙ a1x ∙ a0

Если а3 > 0 тогда возможны следующие случаи:

1)Существует 3 вещественных корня

2)Существует 1 вещественный корень

3) Существует 2 вещественных корня

4) когда дискриминант при нахождении экстремумов функции D=0, то тогда тоже возможен один корень

Для случая а3 < 0 функция будет вести себя также, но сначала она будет убывать.

Для реализации программы и составления блок-схемы необходимо учесть все эти случаи и исследовать функцию: найти экстремумы функции, чтобы понять как она себя ведет на определенных интервалах и задать начальное приближение. Это необходимо для того чтобы выбрать подходящий метод поиска корней.

2.1.1 Метод секущих (метод хорд)

Если x0, x1- приближенные значения корня уравнения f(x) = 0 , а f(x0)f(x1)<0, то последующие приближения находят по формуле:

x_n+1=x_n-(f(x_n)/(f(x_n)-f(x_n-1)))*(x_n-x_n-1), n=1,2,...

Методом хорд называют также метод, при котором один из концов отрезка[a; b] закреплен, т. е. вычисление приближения корня уравнения f(x) = 0 производят по формулам:

x₀=b, x_n+1=x_n-(f(x_n)/(f(x_n)-f(a)))*(x_n-a)

либо

x₀=a, x_n+1=x_n-(f(x_n)/(f(x_n)-f(b)))*(x_n-b)

При этом предполагается, что корень уравнения находится на отрезке[a; b], а f''(x) сохраняет знак на [a; b].

2.1.2 Метод касательных

Если x0 - начальное приближение корня уравнения f(x) = 0, то последовательные приближения находят по формуле

x_n+1=x_n-(f(x_n)/f'(x_n)), n=0,1,2,...

Если f' и f'' непрерывны и сохраняют определенные знаки на отрезке[a; b] , а f(a)f(b) < 0 , то, исходя из начального приближения x0 принадлежит [a; b] удовлетворяющего условию f(x0)f''(x0)>0, можно вычислить с любой точностью единственный корень уравнения f(x) = 0.

2.2 Итерационные методы поиска корней полинома.

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

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

f(x_n+h) =f(x_n)+hf '(x_n) + (h2/2)f "(x_n) +...

Члены, содержащие h во второй и более высоких степенях, отбрасывают; используется соотношение xn+h=xn+1. Предполагается, что переход от хn к хn+1 приближает значение функции к нулю так, что f(xn+h)=0. Тогда

x_n+1=x_n-f(x_n)/f'(x_n)

Значение х_n+1 соответствует точке, в которой касательная к кривой в точке хn пересекает ось х. Так как кривая f(x) отлична от прямой, то значение функции f(x_n+l) скорее всего не будет в точности равно нулю. Поэтому всю процедуру повторяют, причем вместо х_n используют х_n+1. Счет прекращается по достижении достаточно малого значения f(x_n+l). Совершенно ясно, что быстрота сходимости в большой мере зависит от удачного выбора исходной точки. Если в процессе итераций тангенс f'(х) угла наклонной касательной обращается в нуль, то применение метода усложняется. Можно также показать, что в случае бесконечно большого f "(x) метод также не является достаточно эффективным. Так как условие кратности корней имеет вид f(х)=f'(*)=0, то в этом случае метод Ньютона не обеспечивает сходимость. Отметим, что иногда используют другой способ контроля сходимости, состоящий в сравнении х_n и х_n+1.

Программа на С++ для определения корней полинома 3-ей степени по методу Ньютона

#include <conio.h>;

#include <stdio.h>

#include <iostream.h>;

#include <stdlib.h>;

#include <math.h>;

double New(double,double,double,double,double,double);

void Polinom (double,double,double,double,double,double []); /*Объявления, так же называемое прототипом функции*/

void main()

{

clrscr();

double A,B,C,D,R,y;

double Reh[4];

int t;

printf("Введите полином вида: A*X^4+B*x^3+C*X^2+D*X+R \n A=" );

cin>>A;

printf("\nB=");

cin>>B;

printf("\n C=");

cin>>C;

printf("\n D=");

cin>>D;

printf("\n R=");

cin>>R;

printf("\n Введите точность, вида 0.00..01 \n y=");

cin>>y;

if(A==0 && B!=0) /*Т.е. если ввели полином 3-й степени*/

{

Polinom(B,C,D,R,y,Reh);

t=Reh[0];

switch(t)

{

case 1: printf("\n Kolichestvo korney =3,\n Koren 1=%f \n Koren 2= %f \n Koren 3= %f",Reh[1],Reh[2],Reh[3]);

break;

case 2:case 4:case 7: printf("\n Kolichestvo korney =1,\n Koren 1=%f \n", Reh[1]);

break;

case 3:case 5: printf("\n Kolichestvo korney =2,\n Koren 1=%f \n Koren 2= %f \n ",Reh[1],Reh[2]);

break;

case 8: printf("\n Korney net");

break;

}

}

else if(a!=0)

{

/*Расчет полинома 4-й степени, где так же будет использоваться подпрограмма Polinom*/

}

else printf("Введены некоректные данные");

}

/*Подпрограмма для решения полинома 3-й степени*/

void Polinom (double A1,double B1,double C1,double D1,double y,double Kor[4])

{

double X1,X2,D,T,F1,F2;

int t;

D=4*B1*B1-4*3*A1*C1; /*Считаем дискрименант полинома 3-й степени*/

if(D>=0)

{

X1=(-2*B1+sqrt(D))/(6*A1);

X2=(-2*B1-sqrt(D))/(6*A1);

if(X1>X2) /*Сортировка Х1 и Х2 по возростанию*/

{ T=X1;

X1=X2;

X2=T;

}

F1=A1*X1*X1*X1+B1*X1*X1+C1*X1+D1;

F2=A1*X2*X2*X2+B1*X2*X2+C1*X2+D1;

Kor[0]=0;

/*Нахдим количество корней и вид графика функции,

при этом расматриваются все возможные типы построения графика (для наших целей они разбиваются на 8 групп)

и № типа заносится в переменную Kor[0] */

if(A1>0)

if(F1>0)

if(F2<0)

Kor[0]=1;

else if(F2>0)

Kor[0]=2;

else Kor[0]=3;

else if (F1<0)

Kor[0]=4;

else Kor[0]=5;

else

if(F1<0)

if(F2>0)

Kor[0]=1;

else if(F2<0)

Kor[0]=2;

else Kor[0]=3;

else if (F1>0)

Kor[0]=4;

else Kor[0]=5;

if (D==0)

if((F1>0 && A1<0)||(F1<0 && A1>0)) Kor[0]=2;

else

if(F1==0) Kor[0]=7;

else Kor[0]=8;

t=Kor[0]; /*Необходимо т.к. функция switch работает только с типом int*/

switch(t) {

case 1: Kor[1]=New(A1,B1,C1,D1,X1,-y);

X1=(X1+X2)/2;

Kor[2]=New(A1,B1,C1,D1,X1,y);

Kor[3]=New(A1,B1,C1,D1,X2,y);

break;

case 2: Kor[1]=New(A1,B1,C1,D1,X1,-y);

break;

case 3: Kor[1]=New(A1,B1,C1,D1,X1,-y);

Kor[2]=X2;

break;

case 4: Kor[1]=New(A1,B1,C1,D1,X2,y);

break;

case 5: Kor[1]=X1;

Kor[2]=New(A1,B1,C1,D1,X2,y);

break;

case 7: Kor[1]=X1;

break;

}

}

else printf("D<0, производная полинома имеет комплексные корни");

}

/*Функция New предназначена для нахождения корня с необходимой точностью*/

double New(double A1,double B1,double C1,double D1,double G,double e)

{

double P[1000];

/*резервирование максимального количества ячеек памяти,

так как нельзя предугадать количество шагов для нахождения

корня с необходимой точностью*/

double prov;

int i,S;

P[0]=G+e;

S=0;

i=0;

while(S==0)

{

P[i+1]=P[i]-(A1*P[i]*P[i]*P[i]+B1*P[i]*P[i]+C1*P[i]+D1)/(3*A1*P[i]*P[i]+2*B1*P[i]+C1);

prov=fabs(fabs(P[i])-fabs(P[i+1]));

if(prov>fabs(e))

i=i+1;

else S=1;

}

return P[i+1];

}