Как использовать функцию scipy.brentq, чтобы сделать ее эквивалентной коду pascal?

Question

Я пытаюсь заняться археологией программного обеспечения. Попытка заменить старый файл pascal современным скриптом в python.

В какой-то момент скрипт pascal использует процедуру, которая выглядит так:

FUNCTION zbrent(x1,x2,tol: real;fx:Func): real;
(* in the range x1 to x2, ZBRENT searches for a zero point of fx(x:real):real
  until a accuracy of tol is reached. fx must change sign in [x1,x2].
  fx must be defined with the far call option on $F+, and not be nested *)
LABEL 99;
CONST
   itmax=100;
   eps=3.0e-8;
VAR
   a,b,c,d,e: real;
   min1,min2,min: real;
   fa,fb,fc,p,q,r: real;
   s,tol1,xm: real;
   iter: integer;
BEGIN
   a := x1;
   b := x2;
   fa := fx(a);
   fb := fx(b);
   IF (fb*fa > 0.0) THEN BEGIN
      GotoXY(EelX,EelY+8);
      writeln('pause in routine ZBRENT');
      write('Root must be bracketed');
      Readln;
      GotoXY(EelX,EelY+8);
      Writeln('                       ');
      Write('                      ');
      Goto 99;
    END;
   fc := fb;
   FOR iter := 1 to itmax DO BEGIN
      IF (fb/abs(fb)*fc > 0.0) THEN BEGIN
         c := a;
         fc := fa;
         d := b-a;
         e := d
      END;
      IF (abs(fc) < abs(fb)) THEN BEGIN
         a := b;
         b := c;
         c := a;
         fa := fb;
         fb := fc;
         fc := fa
      END;
      tol1 := 2.0*eps*abs(b)+0.5*tol;
      xm := 0.5*(c-b);
      IF ((abs(xm) <= tol1) OR (fb = 0.0)) THEN BEGIN
         zbrent := b; GOTO 99 END;
      IF ((abs(e) >= tol1) AND (abs(fa) > abs(fb))) THEN BEGIN
         s := fb/fa;
         IF (a = c)  THEN BEGIN
            p := 2.0*xm*s;
            q := 1.0-s
         END ELSE BEGIN
            q := fa/fc;
            r := fb/fc;
            p := s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
            q := (q-1.0)*(r-1.0)*(s-1.0)
         END;
         IF (p > 0.0) THEN  q := -q;
         p := abs(p);
         min1 := 3.0*xm*q-abs(tol1*q);
         min2 := abs(e*q);
         IF (min1 < min2) THEN min := min1 ELSE min := min2;
         IF (2.0*p < min) THEN BEGIN
            e := d;
            d := p/q
         END ELSE BEGIN
            d := xm;
            e := d
         END
      END ELSE BEGIN
         d := xm;
         e := d
      END;
      a := b;
      fa := fb;
      IF (abs(d) > tol1) THEN BEGIN
         b := b+d
      END ELSE BEGIN
         IF (xm > 0) THEN BEGIN
            b := b+abs(tol1)
         END ELSE BEGIN
            b := b-abs(tol1)
         END
      END;
      fb := fx(b)
   END;
   writeln('pause in routine ZBRENT');
   writeln('maximum number of iterations exceeded'); readln;
   zbrent := b;
99:   END;

Он описан в этой книге, на странице 285. Это метод Ван Вейнгаардена-Деккера-Брента.

Я хочу заменить его одной строкой в ​​python, желательно используя scipy . Я вижу, что существует метод scipy.optimize.brentq , но есть одно существенное отличие:

• Pascal zbrent использует только один вход допуска (tol), а python имеет rtol и xtol. Я не понимаю, что они означают.

Что с этим делать? Что я должен указать в качестве xtol и rtol в моей программе python, чтобы сделать ее эквивалентной программе pascal? Я ничего не знаю о численных вычислениях. Я напуган. Я просто ученый-материаловед!

Answer 1

Допуск tol в коде Pascal применяется однозначно к длине оставшегося интервала брекетинга как абсолютный допуск. Его функция соответствует параметру xtol в scipy. Вы можете игнорировать rtol и оставить значение по умолчанию. В стандартных ситуациях это не должно иметь значения.

В общем, для больших корней и, соответственно, больших значений для конечных точек интервала брекетинга достижимая или желаемая точность будет иной, чем если бы эти значения были меньше. Это можно контролировать с помощью соответствующим образом масштабированных значений для xtol, если начальный интервал достаточно узкий. Однако, если начальный интервал по какой-то причине равен [1e-6,1e6], тогда может быть целесообразно контролировать точность вывода в основном через относительный допуск.

Из документации brentq

xtol : number, optional
rtol : number, optional
    The computed root ``x0`` will satisfy 
        ``np.allclose(x, x0, atol=xtol, rtol=rtol)``, 
    where ``x`` is the exact root. 

    For nice functions, Brent's method will often satisfy the 
    above condition with ``xtol/2`` and ``rtol/2``. [Brent1973]_

и от np.allclose(a, b, rtol=.., atol=..)

The tolerance values are positive, typically very small numbers.  The
relative difference (`rtol` * abs(`b`)) and the absolute difference
`atol` are added together to compare against the absolute difference
between `a` and `b`.