home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / RTSAFE.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  1.5 KB  |  61 lines
  1. FUNCTION rtsafe(x1,x2,xacc: real): real;
  2. (* Programs using routine RTSAFE must externally define procedure
  3. funcd(x,f,df:real) which returns the function value f and its
  4. derivative df at the point x. *)
  5. LABEL 99;
  6. CONST
  7.    maxit=100;
  8. VAR
  9.    df,dx,dxold,f,fh,fl: real;
  10.    swap,temp,xh,xl,rts: real;
  11.    j: integer;
  12. BEGIN
  13.    funcd(x1,fl,df);
  14.    funcd(x2,fh,df);
  15.    IF (fl*fh >= 0.0) THEN BEGIN
  16.       writeln('pause in routine RTSAFE');
  17.       writeln('root must be bracketed'); readln
  18.    END;
  19.    IF (fl < 0.0) THEN BEGIN
  20.       xl := x1;
  21.       xh := x2
  22.    END ELSE BEGIN
  23.       xh := x1;
  24.       xl := x2;
  25.       swap := fl;
  26.       fl := fh;
  27.       fh := swap
  28.    END;
  29.    rts := 0.5*(x1+x2);
  30.    dxold := abs(x2-x1);
  31.    dx := dxold;
  32.    funcd(rts,f,df);
  33.    FOR j := 1 TO maxit DO BEGIN
  34.       IF((((rts-xh)*df-f)*((rts-xl)*df-f) >= 0.0)
  35.       OR (abs(2.0*f) > abs(dxold*df))) THEN BEGIN
  36.          dxold := dx;
  37.          dx := 0.5*(xh-xl);
  38.          rts := xl+dx;
  39.          IF (xl = rts) THEN GOTO 99 END
  40.       ELSE BEGIN
  41.          dxold := dx;
  42.          dx := f/df;
  43.          temp := rts;
  44.          rts := rts-dx;
  45.          IF (temp = rts) THEN GOTO 99
  46.       END;
  47.       IF (abs(dx) < xacc) THEN GOTO 99;
  48.       funcd(rts,f,df);
  49.       IF (f < 0.0)  THEN BEGIN
  50.          xl := rts;
  51.          fl := f
  52.       END ELSE BEGIN
  53.          xh := rts;
  54.          fh := f
  55.       END
  56.    END;
  57.    writeln('pause in RTSAFE');
  58.    writeln('maximum number of iterations exceeded'); readln;
  59. 99:   rtsafe := rts
  60. END;
  61.