Данная статья посвящена нахождению корней уравнения

        F(x)=0,    (1)

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

function Metod1(f: TFunc; x0: double): double;

В данном методе для удобства вычислений переходят от исходного уравнения (1) к уравнению

        x=f(x).                                         (2)

Данный переход можно осуществить множеством способов, например, прибавить к обеим частям (1) x.
Суть метода последовательных приближений состоит в том, что начальное приближение x[0] подставляется в правую часть формулы (2) и вычисляется значение x[1]. Затем полученное x[1] снова подставляется в правую часть формулы (2) и вычисляется x[2], потом x[3] и т.д. Рабочая формула метода последовательных приближений имеет вид

        x[n]=f(x[n-1]).                                 (3)

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

        | x[n]-x[n-1] | < Eps.                          (4)

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

        | f'(x[n]) | < 1                                (5)

для всех значений x[n].

function Metod2(f: TFunc; x0: double): double;

Формула данного итерационного метода имеет вид

        x[n+1] = x[n] + A*( f(x[n]) - x[n] ),           (6)

где A определяется по формулам

        A = 1/(1 - f'(s))                               (7)
        f'(s) = ( f(x[n]) - x[n] )/( x[n] - x[n-1] )    (8)

при этом на первом шаге x[1] определяется простым методом последовательных приближений.

function Metod3(f: TFunc; x0: double): double;

Небольшая дальнейшая модификация усовершенствованного метода последовательных приближений приводит к одному из наиболее известных численных методов решения уравнений - методу Ньютона-Рафсона. Формула метода для f(x), подчиняющегося соотношению (2), имеет вид

        x[n+1] = (f(x[n])-x[n]*f'(x[n]))/(1-f'(x[n])),  (9)

при этом сходимость метода обеспечивается, если:

1. x[0] выбрано достаточно близко к решению x = f(x);
2. Вторая производная от f(x) не становится слишком большой;
3. Производная f'(x) не слишком близка к 1.

Формула Ньютона-Рафсона для F(x), подчиняющегося соотношению (1), имеет вид

        x[n+1] = x[n] - F(x[n])/F'(x[n]),               (10)

при этом условия сходимости принимают вид:

1. x[0] выбрано достаточно близко к решению x = F(x);
2. Вторая производная f''(x) не становится слишком большой;
3. Производная f'(x) не слишком близка к 0.

Применим метод Ньютона-Рафсона согласно формуле (10), при этом вычисление F(x) будем осуществлять по правилу Горнера

        b[m] = a[m]
        b[j] = a[j] + x[n]*b[j+1], j=m-1, ..., 0        (11)

Таким образом находим F(x) = b[0]. F'(x) представляет собой многочлен степени m-1. Воспользовавшись для его вычисления теми же рекуррентными формулами, имеем:

        c[m] = b[m]
        c[j] = b[j] + x[n]*c[j+1], j=m-1, ..., 1        (12)

и соответственно F'(x) = c[1]. Подставляя найденные значения F(x) и F'(x) в формулу (10) для метода Ньютона-Рафсона, получаем

        x[n+1] = x[n] - b[0]/c[1],                      (13)

где b[o] и c[1] вычислены по формулам (11) и (12). Тихо и незаметно мы пришли к методу Бриге-Виетта ;)

{
  DonNTU. e-moe aka Labinskiy Nikolay (c) 2005
  MailTo: e-moe@ukr.net and visit www.sources.ru
}
 
{$N+,G+,X+}
const
  Eps = 0.00001;        { Точность вычислений }
  anMax = 3;            { Степень полинома }
  { Коэффициенты полинома F(x) = a[0]+a[1]*x+a[2]*x^2+ ... +a[n]*x^n }
  An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13);
                                  { a[0]   a[1]  a[2]  a[3] }
type
  { Описание функции для передачи ее в кач-ве параметра функциям }
  TFunc = function(x: double; px: byte; d: shortint): double;
var
  x0,x1,x2: double; { Переменные для вычисления корней }
 
function Func(x: double; px: byte; d: shortint): double; far;
(*
  Функция для вычисления полинома в точке 'x'
  x  - в этой точке будет происходить расчет
  px - если 1, то результат будет = F(x)+x
  d  - если 0, то F(x), 1 - F'(x), -1 - F(x)/F'(x)
*)
var
  bn,                   { Это будет F(x) }
  cn: double;           { Это будет F'(x) }
  j: byte;              { Счетчик цикла ;) }
begin
  if px=1 then          { Если px=1, то добавляем к рез-ту x }
    An[1]:=An[1]+1;
{ Далее идет расчет по правилу Горнера  }
  bn:=An[anMax];
  cn:=bn;
  for j:= anMax-1 downto 1 do
    begin
      bn:=An[j]+x*bn;
      cn:=bn+x*cn;
    end;
  bn:=An[0]+x*bn;
{ В зависимости от d возвращаем результат }
  case d of
    -1: Func:=bn/cn;
     0: Func:=bn;
     1: Func:=cn;
  else
    Func:=bn;
  end;
{ Забираем обратно свой x }
  if px=1 then
    An[1]:=An[1]-1;
end;
 
function Sign(x: double): shortint;
{ Функция определения знака }
begin
  if x < - Eps then
    Sign:=-1
  else
    if x > Eps then
      Sign:=1
    else
      Sign:=0;
end;
 
function Get_X0(f: TFunc): double;
(*
  Поиск начального приближения:
  Идет начальный просмотр функции и определение двух точек,
  между которыми функция меняят знак (проходит через 0).
  Где-то между этими точками лежит корень.
*)
const
  LookFrom = -10.0;     { Поиск начинается с этой точки }
  step = 1;             { Шаг, с которым производится поиск }
begin
  x0:=LookFrom;
  x1:=x0+eps;
  while Sign(f(x1,0,0)) = Sign(f(x0,0,0)) do
    begin
      x0:=x1;
      x1:=x1+step;
    end;
  Get_X0:= (x1+x0)/2;
end;
 
function Metod1(f: TFunc; x0: double): double;
var
  cnt: word;
function TestFunc(x: double): boolean;
begin
  if abs(f(x,1,1)) >= 1 then
    begin
      writeln('  Ошибка: Невозможно расчитать результат!');
      writeln('  Модуль производной в точке ',x:0:6,' больше 1',#10);
      TestFunc:= false;
    end
  else
    TestFunc:= true;
end;
begin
  Writeln('Метод последовательных приближений:');
  if TestFunc(x0) then
    x1:=f(x0,1,0)
  else
    exit;
{  x2:=f(x1,0)+x1;}
  cnt:=1;
  while abs(x1-x0) > eps do
    begin
      x0:=x1;
      if TestFunc(x0) then
        x1:=f(x0,1,0)
      else
        exit;
      inc(cnt)
    end;
  Metod1:=x1;
  Writeln('x=', x1:0:6);
  Writeln('Метод сошелся на ',cnt,' шаге.',#10);
end;
 
function Metod2(f: TFunc; x0: double): double;
var
  alph: double;
  cnt: word;
begin
  Writeln('Усовершенствованный метод последовательных приближений:');
  x1:=f(x0,1,0);
  alph:=1/(1-f(x1,0,0)/(x1-x0));
  x2:=x1+alph*f(x1,0,0);
  cnt:=2;
  while abs(x2-x1) > eps do
    begin
      x0:=x1;
      x1:=x2;
      alph:=1/(1-f(x1,0,0)/(x1-x0));
      x2:=x1+alph*f(x1,0,0);
      inc(cnt)
    end;
  Metod2:=x2;
  Writeln('x=', x2:0:6);
  Writeln('Метод сошелся на ',cnt,' шаге.',#10);
end;
 
function Metod3(f: TFunc; x0: double): double;
var
  cnt: word;
begin
  Writeln('Метод Бирге-Виетта:');
  x1:=x0-f(x0,0,-1);
  cnt:=1;
  while abs(x1-x0) > Eps do
    begin
      x0:=x1;
      x1:=x0-f(x0,0,-1);
      inc(cnt);
    end;
  Metod3:=x1;
  Writeln('x=', x1:0:6);
  Writeln('Метод сошелся на ',cnt,' шаге.',#10);
end;
 
begin
  x0:=Get_X0(Func);
  Writeln;
  Writeln('F(x)= -5.372 + 1.2493*x + 0.559*x^2 - 0.13*x^3');
  Writeln('Eps = ',eps:0:5,'  x0= ',x0:0:6,#10);
 
  Metod1(Func,x0);
  Metod2(Func,x0);
  Metod3(Func,x0);
 
  Write('Press Enter to exit');
  Readln;
end.

Автор: e-moe aka Labinskiy Nikolay © 2005, DonNTU