====== Численное решение уравнений итерационными методами ======

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

==== Метод последовательных приближений ====
<code pascal>function Metod1(f: TFunc; x0: double): double;</code>
В данном методе для удобства вычислений переходят от исходного уравнения (1) к уравнению
<code>        x=f(x).                                         (2)</code>
Данный переход можно осуществить множеством способов, например, прибавить к обеим частям (1) //x//.\\ Суть метода последовательных приближений состоит в том, что начальное приближение //x[0]// подставляется в правую часть
формулы (2) и вычисляется значение //x[1]//. Затем полученное //x[1]// снова подставляется в правую часть формулы (2)
и вычисляется //x[2]//, потом //x[3]// и т.д. Рабочая формула метода последовательных приближений имеет вид<code>        x[n]=f(x[n-1]).                                 (3)</code>
Вычисления продолжаются до тех пор, пока не будет достигнута заданная точность //Eps//, т.е.
<code>        | x[n]-x[n-1] | < Eps.                          (4)</code>
Основной проблемой при работе с итерационными методами является **проблема сходимости**. Достаточным условием сходимости метода последовательных приближений является выполнение условия
<code>        | f'(x[n]) | < 1                                (5)</code>
для всех значений x[n].

==== Усовершенствованный метод последовательных приближений ====
<code pascal>function Metod2(f: TFunc; x0: double): double;</code>
Формула данного итерационного метода имеет вид
<code>        x[n+1] = x[n] + A*( f(x[n]) - x[n] ),           (6)</code>
где A определяется по формулам
<code>        A = 1/(1 - f'(s))                               (7)
        f'(s) = ( f(x[n]) - x[n] )/( x[n] - x[n-1] )    (8)</code>
при этом на первом шаге //x[1]// определяется простым методом последовательных приближений.

==== Метод Ньютона-Рафсона, Бриге-Виетта ====
<code pascal>function Metod3(f: TFunc; x0: double): double;</code>
Небольшая дальнейшая модификация усовершенствованного метода последовательных приближений приводит к одному из
наиболее известных численных методов решения уравнений - //методу Ньютона-Рафсона//. Формула метода для //f(x)//, подчиняющегося соотношению (2), имеет вид
<code>        x[n+1] = (f(x[n])-x[n]*f'(x[n]))/(1-f'(x[n])),  (9)</code>
при этом //сходимость метода обеспечивается//, если:

  - //x[0]// выбрано достаточно близко к решению //x = f(x)//;
  - Вторая производная от //f(x)// не становится слишком большой;
  - Производная //f'(x)// не слишком близка к 1.
Формула Ньютона-Рафсона для //F(x)//, подчиняющегося соотношению (1), имеет вид
<code>        x[n+1] = x[n] - F(x[n])/F'(x[n]),               (10)</code>

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

Применим метод Ньютона-Рафсона согласно формуле (10), при этом вычисление //F(x)// будем осуществлять по //правилу Горнера//
== с использованием рекуррентных формул: ==
<code>        b[m] = a[m]
        b[j] = a[j] + x[n]*b[j+1], j=m-1, ..., 0        (11)</code>
Таким образом находим //F(x) = b[0]//.
F'(x) представляет собой многочлен степени //m-1//. Воспользовавшись для его вычисления теми же рекуррентными формулами, имеем:
<code>        c[m] = b[m]
        c[j] = b[j] + x[n]*c[j+1], j=m-1, ..., 1        (12)</code>
и соответственно //F'(x) = c[1]//.
Подставляя найденные значения //F(x)// и //F'(x)// в формулу (10) для метода Ньютона-Рафсона, получаем<code>        x[n+1] = x[n] - b[0]/c[1],                      (13)</code>
где //b[o]// и //c[1]// вычислены по формулам (11) и (12).
Тихо и незаметно мы пришли к //методу Бриге-Виетта// ;)

<code pascal>{
  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.</code>

Автор:
e-moe aka Labinskiy Nikolay (c) 2005, DonNTU