unit Equ2;
// Методы поиска корней нелинейного уравнения

interface
uses UMatrix;
// AX0: динамический массив начальных значений
//          должен заполняться внешней процедурой
// DAbsX: заданная абсолютна допустимая погрешность X
// DRelX: заданная относительная допустимая погрешность X
// DF: заданная допустимая невязка
// NIterMax: число обращений к функции, после которого
//           пользователю делается запрос о продолжении поиска
// AX: динамический массив текущих значений переменных
// AFDelta: динамический массив значений невязок

var AX, AX0, AFDelta{, AXOld}: TVector;
    DAbsX: real = 0.001;
    DRelX: real = 0.001;
    DF: real = 0;
    NIter:longint;    // Число обращений к функции
    NIterMax: longint = 10000;

// Func: переменная - функция
// X - заданный вектор переменных, N - его размер,
// Y - вектор правых частей уравнений, рассчитываемый функцией;
// Функция должна также увеличивать на 1 переменную NIter;
// Если в процессе расчета возникли ошибки,
// функция должна возвращать false
var Func: function(X: TVector; var Y: TVector; N: word): boolean;

// Расчет величины допустимой погрешности по x[Ind]
function Deltai(Ind: integer): real;

// Метод простых итераций
function Iterat(N: integer): boolean;

// Модифицированный метод Ньютона
function Nevton(N: integer): boolean;

implementation
uses Forms, Windows, Sysutils, Math;

var AXOld: TVector; // Предыдущие значения переменных
    First: boolean; // Флаг первой итерации

function Deltai(Ind: integer): real;
// Расчет величины допустимой погрешности по x[Ind]
begin
 Result := Max(DAbsX, Abs(AX[Ind] * DRelX));
end;

function Iterat(N: integer): boolean;
var i: word;
    SumFDelta, DX: real;
begin
 Result := false;
 AX := CreateVector(N);
 AFDelta := CreateVector(N);
 AXOld := CreateVector(N);
 AX := Copy(AX0, 0, N);
 NIter := 0;
 First := true;
 repeat
  AXOld := Copy(AX, 0, N);
  if not Func(AX, AX, N ) then exit;
// Вектор AFDelta приращений переменных и его модуль SumFDelta:
  AFDelta := SubVector(AX, AXOld);
  SumFDelta := ModVector(AFDelta);
  // Если хотя бы по одной координате
  // условия окончания не выполняются,
  // расчет продолжается
  Result := true;
  for i := 0 to N - 1 do
   if (Deltai(i) <= Abs(AFDelta[i]))
    then begin
     Result := false;
     break;
    end;
// Окончание
  if Result then exit;
// Проверка сходимости
  if not First
   then if (SumFDelta > DX)
    then begin
     Application.MessageBox(
                   'Метод расходится!' + #13 +
                   'Смените исходную точку или метод',
                   'Метод простых итераций расходится', MB_OK + MB_ICONSTOP);
     exit;
    end;
  if (NIter mod NIterMax) = 0
   then if Application.MessageBox(
                   PChar('Сделано ' + IntToStr(NIter) + ' итераций' +
                   #13'Продолжать?'),
                   'Подтвердите продолжене расчета',
                   MB_YESNOCANCEL + MB_ICONQUESTION)
                   <> IDYES
    then exit;
  DX := SumFDelta;
  First := false;
 until false;
end;

function Nevton(N: integer): boolean;
var Delta: TVector;    // Поправка Ньютона
    Fk: TVector;       // Вектор невязок при приращении xi
    Jac: TMatrix;      // Матрица Якоби
    FMaxK, FOld: real;// Текщая и предыдущая невязки
    t: real;
    i: integer;

function Jacob: boolean;
// Формирование матрицы Якоби
var i, k: integer;
    D: real;
begin
 Result := true;
 AX := Copy(AXOld, 0, N);
 for i:=0 to N-1 do begin
  D := Deltai(i);
  AXOld[i] := AXOld[i] + D;
  if not Func(AXOld, Fk, N )
   then begin
    Result := false;
    exit;
   end;
  for k:=0 to N-1 do
   Jac[k,i] := (Fk[k] - AFDelta[k]) / D;
  AXOld[i] := AX[i];
 end;
end;

begin
 Result := false;
 AX := CreateVector(N);
 AX := Copy(AX0, 0, N);        // Начальная точка
 AXOld := CreateVector(N);
 AFDelta := CreateVector(N);
 Fk := CreateVector(N);
 Jac := CreateMatrix(N, N);
 NIter := 0;
 if not Func(AX, AFDelta, N)   // Невязки
  then exit;
 FMaxK := ModVector(AFDelta);  // Норма невязки
 First := true;
 repeat
  t := 1;
  AXOld := Copy(AX, 0, N);  // Запоминание начальной точки
  FOld := FMaxK;
  if not Jacob then exit; // Формирование матрицы Якоби
  Delta := LUMaxSystem(Jac, AFDelta);
  AX := SubVector(AXOld, CMultVector(t,Delta));
  if not Func(AX, AFDelta, N)   // Невязки
   then exit;
  FMaxK := ModVector(AFDelta);  // Норма невязки
// Проверка сходимости
  if not First then
   while (FOld < FMaxK) do
   begin
    t := t / 2;
    if (t < 0.001)
     then begin
      Application.MessageBox(
                   'Метод расходится!' + #13 +
                   'Смените исходную точку или метод',
                   'Метод Ньютона расходится', MB_OK + MB_ICONSTOP);
      Result := false;
      exit;
     end;
    AX := SubVector(AXOld, CMultVector(t,Delta));
    if not Func(AX, AFDelta, N)  // Невязки
     then exit;
    FMaxK := ModVector(AFDelta);  // Норма невязки
   end;
// Проверка критериев окончания
  if FMaxK <= DF
   then Result := true
  else begin
   Result := true;
   for i:= 0 to N-1 do
    if (Deltai(i) <= Abs(Delta[i]))
     then begin
      Result := false;
      break;
     end;
  end;
// Окончание
  if Result then exit;
  if (NIter mod NIterMax) = 0
   then if Application.MessageBox(
                   PChar('Сделано ' + IntToStr(NIter) + ' итераций' +
                   #13'Продолжать?'),
                   'Подтвердите продолжене расчета',
                   MB_YESNOCANCEL + MB_ICONQUESTION)
                   <> IDYES
    then begin
     Result := false;
     exit;
    end;
  First := false;
 until false;
end;

end.
