Реализация некоторых численных методов

8.1 Программирование методов решения нелинейных уравнений в Maxima

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

Пусть дано нелинейное уравнение , и необходимо решить это уравнение, т. е. найти его корень .

Если функция имеет вид многочлена степени , где — коэффициенты многочлена, , то уравнение имеет m корней (основная теорема алгебры).

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

Как известно, не всякое уравнение может быть решено точно. В первую очередь это относится к большинству трансцендентных уравнений.

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

Однако точное решение уравнения не всегда является необходимым. Задачу отыскания корней уравнения можно считать практически решенной, если мы сумеем найти корни уравнения с заданной степенью точности. Для этого используются приближенные (численные) методы решения.

Большинство употребляющихся приближенных методов решения уравнений являются, по существу, способами уточнения корней. Для их применения необходимо знание интервала изоляции , в котором лежит уточняемый корень уравнения.

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

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

Процесс определения корней алгебраических и трансцендентных уравнений состоит из двух этапов:

• отделение корней, — т.е. определение интервалов изоляции , внутри которого лежит каждый корень уравнения;
• уточнение корней, — т.е. сужение интервала до величины равной заданной степени точности .

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

• метод половинного деления (метод дихотомии);
• метод простых итераций;
• метод Ньютона (метод касательных);
• модифицированный метод Ньютона (метод секущих);
• метод хорд и др.

8.1.1 Метод половинного деления

Рассмотрим следующую задачу: дано нелинейное уравнение , необходимо найти корень уравнения, принадлежащий интервалу , с заданной точностью .

Для уточнения корня методом половинного деления последовательно осуществляем следующие операции:

• вычисляем значение функции в точках и ;
• если , то корень находится в левой половине интервала , поэтому отбрасываем правую половину интервала и принимаем ;
• если условие не выполняется, то корень находится в правой половине интервала , поэтому отбрасываем левую половину интервала за счёт присваивания .

В обоих случаях новый интервал в 2 раза меньший предыдущего.

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

Реализация метода половинного деления в виде функции Maxima представлена в следующем примере:

• собственно функция , в которую передаётся выражение , определяющее уравнение, которое необходимо решить

  (%i1)	bisect(f,sp,eps):=block([a,b],a:sp[1],b:sp[2],p:0,
  while abs(b-a)>eps do (
  	p:p+1, c:(a+b)/2,
  	fa:float(subst(a,x,f)), fc:float(subst(c,x,f)),
  	if fa*fc<0 then b:c else a:c
  ),
  float(c));

  

• последовательность команд, организующая обращение к bisect и результаты вычислений

  (%i2)	f:exp(-x)-x$
  a:-1$ b:2$ eps:0.000001$ xrez:bisect(f,[-2,2],eps)$
  print("Решение ",xrez," Невязка ",subst(xrez,x,f))$

  Решение 0.56714344024658 Невязка - 2.34815726529724610^-7

  

В представленном примере решается уравнение . Поиск корня осуществляется на отрезке [-2, 2] с точностью 0.000001.

Следует отметить особенность программирования для Maxima, заключающуюся в том, что решаемое уравнение задаётся в виде математического выражения (т.е. фактически текстовой строки). Числовое значение невязки уравнения вычисляется при помощи функции subst, посредством которой выполняется подстановка значений x = a или в заданное выражение. Вычисление невязки после решения осуществляется также путём подстановки результата решения в выражение .