unit Newton;

{*******************************************************
 * knihovna pro aproximaci korenu rovnic pomoci        *
 * Newtonovy metody                                    *
 *                                                     *
 * Vytvoril: Adam Husar, xhusar01@stud.fit.vutbr.cz    *
 * Datum posledni upravy: 9. 4. 2003                   *
 *******************************************************}


interface

uses
  FEvalLib, DerivLib;

const
  MAX_APROX_COUNT = 500;

type
  TArrResult = array[0..MAX_APROX_COUNT] of extended;

  TNewton = class
    public
      constructor Create;
      destructor  Close;

      function Aproxime(strExpr: string; a, b: extended; Eps: extended;
                 var arrResult:TArrResult; var Count: integer; var isOK: boolean)
                 : extended;

      function GetFunc: TFunctionEval;
    private
      function CheckConditions(a, b: extended): boolean;

    public
      Func: TFunctionEval;
    private
      Derive: TDerive;
  end;

implementation
uses
  Errors;


constructor TNewton.Create;
begin
  Func:= TFunctionEval.Create;
  Derive:= TDerive.Create(0);
end;

destructor TNewton.Close;
begin
  Func.Close;
  Derive.Close;
end;

function TNewton.GetFunc: TFunctionEval;
begin
  GetFunc:= Func;
end;


function TNewton.CheckConditions(a, b: extended): boolean;
{fce vraci true pokud jsou pocattecni podminky pro krajni body splneny,
  jinak false}
var
  bRes: boolean;
  auxExt: extended;
  da1, da2: integer;
begin
  if a > b then begin
    auxExt:= a;      {vymena}
    a:= b;
    b:= auxExt;
  end;

{maji opacne znamenko?}
  bRes:= Func.EvalExpr(a) * Func.EvalExpr(b) < 0;
  if not bRes then begin
    if GetLastErrorCode = ERR_OK then
      SetLastErrorCode(ERR_N_POINT_SGN_EQ);
    CheckConditions:= bRes;
    exit;
  end;

{ma 1. der. stejne zn.?}
  bRes:= Derive.CountFirst(Func, a) * Derive.CountFirst(Func, b) >= 0;
  if not bRes then begin
    if GetLastErrorCode = ERR_OK then
      SetLastErrorCode(ERR_N_POINT_D1_SGN_NEQ);
    CheckConditions:= bRes;
    exit;
  end;

{ma 2.der stejne znamenko?}
  da2:= Derive.GetSecondSign(Func, a, +1);
  da1:= Derive.GetSecondSign(Func, b, -1);
  bRes:= (da1 = da2) or (da1 = 0) or (da2 = 0);
  if not bRes then begin
    if GetLastErrorCode = ERR_OK then
      SetLastErrorCode(ERR_N_POINT_D2_SGN_NEQ);
    CheckConditions:= bRes;
    exit;
  end;

  CheckConditions:= bRes;
end;



function TNewton.Aproxime(strExpr: string; a, b: extended; Eps: extended;
       var arrResult:TArrResult; var Count: integer; var isOK: boolean): extended;
{funkce vraci koren rovnice vyhovujici zadane presnosti, v poli Result
 postupne aproximace s poctem v Count, jestli vse probehlo dobra je po
 ukonceni isOk true}
var
  i, d2a, d2b: integer; {pomocny integer}
  fd: extended;
  bAproxOk: boolean;
begin

  Derive.SetAccuracy(Eps);  {nastaveni presnosti pro derivace}

  i:= Func.ReadExpr(strExpr); {nacteni fcniho vyrazu}
  isOK:= i = 0;
  if not isOk then begin
    Aproxime:= 0;
    exit;
  end;

  isOK:= CheckConditions(a, b); {kontrola pocatecnich podminek}
  if not isOk then begin
    Aproxime:= 0;
    exit;
  end;

{urceni prvniho bodu:}
  d2a:= Derive.GetSecondSign(Func, a, +1);
  d2b:= Derive.GetSecondSign(Func, b, -1);

  if Func.EvalExpr(a) * d2a > 0 then
    arrResult[0]:= a
  else if Func.EvalExpr(b) * d2b > 0 then
    arrResult[0]:= b
  else if (d2a = 0) or (d2b = 0) then
    arrResult[0]:= a
  else begin
    SetLastErrorCode(ERR_N_UNKN_FIRST_POINT);
    isOk:= false;
    Aproxime:= 0;
    exit;
  end;

{vypocet aproximaci}
  i:= 1;
  repeat
    fd:= Derive.CountFirst(Func, arrResult[i-1]);

    if fd <> 0 then
      arrResult[i]:= arrResult[i-1] - Func.EvalExpr(arrResult[i-1]) / fd
         {vypocet aprox}
    else begin
      SetLastErrorCode(ERR_ZERO_DIV);
      isOK:= false;
      Aproxime:= 0;
      exit;
    end;

    bAproxOk := abs(arrResult[i] - arrResult[i-1]) < Eps;
      {mame uz dostatecnou presn.?}

    i:= i+1;
    if not bAproxOk and  (i >= MAX_APROX_COUNT) then begin
      SetLastErrorCode(ERR_N_BIGGER_APROX_COUNT);
      isOk:= false;
      Aproxime:= arrResult[i-1];
      exit;
    end;

  until bAproxOk or (GetLastErrorCode <> ERR_OK);

  Count:= i-1;
  Aproxime:= arrResult[Count];
end;

end.