Опубликован: 10.12.2015 | Уровень: для всех | Доступ: платный
Лекция 12:

Численное решение алгебраических уравнений

< Лекция 11 || Лекция 12 || Лекция 13 >

Простые линейные или квадратные уравнения легко решить программами для расчета по формулам. Хотя алгебраические уравнения третьей и четвертой степени еще можно решить аналитически, соответствующие формулы настолько сложны, что намного предпочтительнее бывает воспользоваться численными методами.

Принципы численного решения алгебраических уравнений подробно описаны в соответствующей литературе (К. Эберт, Х.Эдерер Компьютеры. Применение в химии. М.: Мир. 1988. с. 109-139; А.А. Самарский, А.В. Гулин Численные методы. М.: Наука. 1989. с. 48-59; Р.С. Гулин, Б.В. Овчинский Элементы численного анализа и математической обработки результатов опыта. М.: Наука. 1970. с. 24-45).

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

Методом Ньютона (метод касательных).

Для наглядности рассмотрим рисунок.


Пусть при поиске корня получено значение xi, (это значение можно считать начальным приближением). Сначала рассчитывают значение функции при x=xi. Затем проводят касательную к графику функции f(x) в точке с координатами xi, f(xi) и находят точку пересечения этой касательной с осью 0x. Точка пересечения рассматривается теперь как новое приближенное значение корня, и описанная процедура повторяется с новым значением х. Как видно из рисунка, точки пересечения касательных с осью х быстро сходятся к точке пересечения графика функции f(x) с осью 0x.

Касательную в точке х можно построить с помощью производной df(x)/dx. Из рассмотрения треугольника, образованного f(xi), х, и xi+1, следует

df(xi)/dx=f'(xi)=f(xi)/(xi-xi+1)

После простых преобразований получается итерационная формула метода Ньютона:

x i+1=xi- f(xi)/f'(i)

Из этой рекуррентной формулы видно, что на каждом итерационном шаге используется как значение функции f(xi), так и значение производной f'(xi). Кроме того, необходимо задать начальное значение x, от которого зависят достоверность результата и время, необходимое для его получения.

Ниже приводится листинг программы. В функции Equation задается уравнение. Производная функции в точке x рассчитывается в функции root и вызывается в функции Derivative, при этом функция Equation вызывается через указатель.

Пользователь задает начальное приближение (t). Далее, переменной y присваивается значение функции в этой начальной точке, а переменной y1 присваивается значение производной в этой точке. Затем вычисляется не слишком ли мало значение производной. Геометрически это равнозначно вопросу, не слишком ли полого расположена касательная. Если значение производной слишком мало, то следует повторный запрос на ввод более приемлемого начального приближения. Затем следует проверка сходимости результата: проверяется совпадают ли в значениях x, полученных в последних двух итерациях, первые 7 знаков. Если совпадают, то управление передается оператору вывода, в противном случае новое значение x присваивается переменной t и итерация повторяется.

//==========================================================
// Name        : equation_newton.cpp
// Author      : Marat
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//==========================================================

#include <iostream>
using namespace std;

#include<math.h>

//Здесь описывается уравнение функции
double Equation (double x)
{
	return (10*pow(x,5)-exp(pow(x,2)));
}

//Здесь производится дифференцирование. Вызов функции
// Equation осуществляется через указатель.
//Определим тип указателя на функцию Equation
//dif - функция дифференцирования
typedef double(*dif)(double);
//Функция расчета производной
double root(dif F, double x)
{
int n=1000000;
double delta=2*x/n;
return ((*F)(x+delta)-(*F)(x))/delta;
}
//Вызов функции дифференцирования
double Derivative(double x)
{	double root(dif F, double x);
	return root(Equation, x);
}

int main() {
	double Equation(double x);
	double Derivative(double x);
	double t;
	alpha: cout<<"\nInput initial value, please ";
	cin>>t;
	gamma: double y=Equation(t);
	double y1=Derivative(t);
	if (fabs(y1)<1.0E-20) goto alpha;
	double t1=t-y/y1;
	if (fabs(t-t1)<fabs(t*1.0E-7)) goto beta;
	t=t1; goto gamma;
	beta: cout<<"\nThe root of equation is "<<(t+t1)/2;
	cout<<"\nFunction value is "<<y<<endl;
	return 0;
}

Результат:

Input initial value, please 1
The root of equation is 0.694941
Function value is 2.23705e-007

Метод половинного деления. Уравнение имеет действительные корни, если найдутся два аргумента x1 и x2, для которых f(x1) и f(x2) имеют разные знаки. Возьмем для примера функцию 0,1x5+x4-3x3+2x2-17x+23=0. Выберем отрезок, на котором функция меняет знак. В нашем случае это интервал [0,2]: f(0) =23 >0, f(2)=-7,8 <0. Следовательно, на отрезке [0,2] имеется хотя бы одно решение. Для нахождения корня уравнения найдем значение функции на середине отрезка: f(1)=6,1 >0. Следовательно, решение находится на интервале [1,2]. С каждым новым шагом этой процедуры отрезок, на котором лежит корень уравнения, уменьшается в два раза. Если отрезок становится достаточно малым, то деление можно прекратить. В рассматриваемом примере корень уравнения на отрезке [0,2] равен 1,34375.


Рассмотренный метод называется методом деления отрезка пополам, поскольку отрезок, содержащий корень уравнения, после очередной итерации делится пополам. Программа предусматривает однократный поиск двух значений аргумента, при которых значения функции имеют разные знаки и итерационную процедуру деления отрезка пополам. Ниже приводится листинг программы.

В приведенном ниже листинге алгоритм решения методом хорд заложен в функцию devide, вызов функции Equation происходит через указатель, поэтому легко существующую программу адаптировать для решения нескольких уравнений, подставляя их имена при вызове функции devide или (chorda, см. ниже) в главной функции.

Метод хорд

Другим часто используемым приближенным методом вычисления корней уравнения является метод хорд (правило секших, правило ложного положения, regula falsi), он очень похож на метод Ньютона. Поиск точки, в которой функция обращается в ноль, начинается с того, что через две точки, лежащие на графике f(x), проводят прямую. Точку пресечения этой прямой с осью ox принимают за новое значение x. Как показано на рисунке, такие построения дают набор лучей:

(f(xi)-0)/(xi-xi+1)=(f(xi-1)-0)/(xi-1-xi+1)

Эта формула отличается от итерационной формулы метода Ньютона тем, что производная в формуле Ньютона заменяется отношением конечных приращений функции и аргумента. При решении уравнений методом хорд требуется вводить два начальных значения независимой переменной, а не одно, как в методе Ньютона.


В приведенном ниже листинге алгоритм решения методом хорд заложен в функцию chorda.

//==========================================================
// Name        : deviding_chorda.cpp
// Author      : Marat
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//==========================================================

#include <iostream>
using namespace std;
#include<math.h>
#include<stdlib.h>

//Здесь задается уравнение
double Equation(double x)
{
	return 0.1*pow(x,5)+pow(x,4)-3*pow(x,3)+2*pow(x,2)-17*x+23;
}
//Определение указателя на функцию Equation
typedef double(*dev)(double);

//Функция решения уравнения методом половинного деления
double devide(dev F, double a,double e)
{
double ep=0.1;//Точность
double h,m,y3;
int n=10;
double y1=(*F)(a);
double y2=(*F)(e);
if (y1*y2<=0) goto alpha;
int i;
gamma: h=(e-a)/n;
for(i=1;i<n-1;i++)
{
	double y9=(*F)(a+i*h);
	if (y9*y1>=0) break;//Выход из цикла
	e=a+i*h;
	goto alpha;
}
n*=10;
if (n<=1000) goto gamma;
cout<<"\nRoots were not founded! ";exit(1);
alpha: //Начало процедуры деления отрезка пополам
y1=(*F)(a);
y2=(*F)(e);
if (fabs(a-e)<fabs(ep)) {return (e+a)/2;}
m=(a+e)/2;
y3=(*F)(m);
if (y3*y1>=0) goto epsilon;
e=m; goto alpha;
epsilon: a=m; goto alpha;
}

//Решение уравнений методом хорд
//(правило пропорциональных частей)
double chorda(dev G, double a,double e)
{
double x;
double x0=a;
double x1=e;
gamma: double y0=(*G)(x0);
double y1=(*G)(x1);
if (fabs(x0-x1)<fabs(x0*1.0E-7))
{
	return (x0+x1)/2;
}
x=x1-(x1-x0)*y1/(y1-y0);
x0=x1;
x1=x;
goto gamma;
}

int main() {
	double a; //Нижний предел
	double e; //Верхний предел
	//Прототип функции половинного деления
	double devide(dev F, double a,double e);
	//Прототип функции решения уравнеия
	//методом хорд
	double chorda(dev G, double a,double e);
	cout<<"\nInput lower level, please, a = ";
	cin>>a;
	cout<<"\nInput upper level, please, e = ";
	cin>>e;
	cout<<"\nRoot made by deviding method is "<<devide(Equation,a,e);
	cout<<"\nRoot made by chords   method is "<<chorda(Equation,a,e)<<endl;
	return 0;
}

Результат:

Input lower level, please, a = 0
Input upper level, please, e = 2
Root made by deviding method is 1.34375
Root made by chords   method is 1.35508

Указатели на функции очень удобны при создании разного рода меню. При этом действия, предлагаемые на выбор пользователю, оформляются в виде функций, вызываемых по указателю. Обращение к функции вызывает выполнение требуемых действий.

В следующем примере рассматривается расчет значения определенного интеграла, выбираемого из меню.

//==========================================================
// Name        : pointer_array_function.cpp
// Author      : Marat
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//==========================================================

#include <iostream>
using namespace std;
#include<math.h>
#include<stdlib.h>

double function_sin(double x)
{
return -cos(x);
}

double function_cos(double x)
{
return sin(x);
}

double function_tg(double x)
{
return -log(fabs(cos(x)));
}

double function_ctg(double x)
{
return log(fabs(sin(x)));
}

double function_pow_x_2(double x)
{
return pow(x,3)/3;
}

typedef double(*menu)(double);
double root(menu F, double A, double B)
{
double s=(*F)(B)-(*F)(A);
return s;
}

int main() {
	double A;
	double B;
	//Нижний предел интегрирования
	cout<<"\nInput lower integrating level, A = ";
	cin>>A;
	//Верхний предел интегрирования
	cout<<"\nInput upper integrating level, B = ";
	cin>>B;

	cout<<"\nChoice the function to be integrated and type number!";
	cout<<"\n1: sin(x) ";
	cout<<"\n2: cos(x) ";
	cout<<"\n3: tg(x)  ";
	cout<<"\n4: ctg(x) ";
	cout<<"\n5: pow(x,2)"<<endl;

	int number;//Номер выбранной функции
	cin>>number;

	double root(menu,double,double);
	double result;

	switch (number)
		{
		case 1:result=root(function_sin,A,B); break;
		case 2:result=root(function_cos,A,B); break;
		case 3:result=root(function_tg, A,B); break;
		case 4:result=root(function_ctg,A,B); break;
		case 5:result=root(function_pow_x_2,A,B); break;
	default: cout<<"\nThis number is illegal!"<<endl;exit(1);
		}
	cout<<"\nThe answer is "<<result<<endl;
	return 0;
}

Результат:

Input lower integrating level, A = 3.5
Input upper integrating level, B = 6.7
Choice the function to be integrated and type number!

1: sin(x)
2: cos(x)
3: tg(x)
4: ctg(x)
5: pow(x,2)

5

The answer is 85.9627

Говоря о вызовах функций через указатели нельзя обойти стороной использование стандартной функции сортировки qsort().

Эта функция имеет следующий прототип:

void qsort(void base, size_t nelem, size_t width, int (*fcmp)(const void *p1, const void *p2));

Прототип находится в заголовочном файле stdlib.h.

Функция qsort() сортирует массив, постоянно вызывая функцию сравнения. Для вызова функции сравнения ее адрес должен заместить указатель fcmp, специфицированный как формальный параметр.

Рассмотрим параметры функции qsort():

base указатель на начало таблицы сортируемых элементов (адрес 0-го элемента массива);

nelem количество элементов в таблице (целая величина, не большая размера массива);

width размер элемента таблицы (целая величина, определяющая в байтах размер одного элемента массива);

fcmp указатель на функцию сравнения, получающую в качестве параметров два указателя p1 и p2 на элементы таблицы и возвращающую в зависимости от результата сравнения целое число: если *pl < *p2, fcmp возвращает целое < 0;

если *pl == *р2, fcmp возвращает 0;

если *pl > *p2, fcmp возвращает целое > 0.

При сравнении символ "меньше чем" (<) означает, что после сортировки левый элемент отношения *p1 должен оказаться в таблице перед правым элементом *р2, т.е. значение *p1 должно иметь меньшее значение индекса в массиве, нежели *р2. Аналогично (но обратно) определяется расположение элементов при выполнении соотношения "больше чем" (>).

В следующей программе функция qsort() упорядочивает массив вещественных чисел сначала по возрастанию, затем – по убыванию.

//==========================================================
// Name        : sort_ascending_descending.cpp
// Author      : Marat
// Version     :
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//==========================================================

#include <iostream>
using namespace std;
#include<stdlib.h>   //Для функции qsort

//Вызов функции упорядочивания:
//(void *)pc-Адрес начала сортируемой таблицы
//n - Количество элементов сортируемой таблицы
//sizeof(pc[0])- Размер одного элемента
//ascending - Имя функции сравнения (указатель) по возрастанию
//descending - Имя функции сравнения (указатель) по убыванию

//Функция сравнения при сортировке по возрастанию
int ascending (const void *a, const void *b)
{double *pa=(double *)a;
 double *pb=(double *)b;
 int c;
 *pa<*pb?c=-1:c=1;
 if(*pa==*pb){c=0;}
 return c;
}

//Функция сравнения при сортировке по убыванию
int descending (const void *a, const void *b)
{double *pa=(double *)a;
 double *pb=(double *)b;
 int c;
 *pa>*pb?c=-1:c=1;
 if(*pa==*pb){c=0;}
 return c;
}

void print(int n, double a[])
{
	for (int i=0;i<n;i++)
{
	cout<<"\na["<<i+1<<"]="<<a[i];
}
cout<<"\n"<<endl;
}

int main() {
	//Зададим вещественный массив
	double pc[]={2.64, 8.16, 0.23, -6.54, 3.25};

	//Определим количество его элементов
	int n=sizeof(pc)/sizeof(pc[0]) ;//Размер таблицы:

	//Выведем на экран неотсортированный массив

	cout<<"\nBefor sorting:";

	void print(int n, double a[]);
	print(n, pc);

	//Отсортируем массив по возрастанию
	qsort((void *)pc, n, sizeof(pc[0]), ascending);

	//Распечатаем массив после сортировки по возрастанию
	cout<<"\nAfter sorting by ascending:";
	print(n, pc);

	//Отсортируем массив по убыванию
	cout<<"\nAfter sorting by descending:";
	qsort((void *)pc, n, sizeof(pc[0]), descending);

	//Расчпечатаем массив после сортировки по убыванию
	print(n, pc);
	return 0;
}

Результат:

Before sorting:

a[1]=2.64
a[2]=8.16
a[3]=0.23
a[4]=-6.54
a[5]=3.25

After sorting by ascending:
a[1]=-6.54
a[2]=0.23
a[3]=2.64
a[4]=3.25
a[5]=8.16

After sorting by descending:
a[1]=8.16
a[2]=3.25
a[3]=2.64
a[4]=0.23
a[5]=-6.54
< Лекция 11 || Лекция 12 || Лекция 13 >
Зося Ковалева
Зося Ковалева

Хочу получить удостоверение. Сколько стоит оплата?

Aleksey Aplaev
Aleksey Aplaev
Россия, Chelybinsk
Александр Сидоров
Александр Сидоров
Россия, Самара