如何使用scipy.brentq函数，使其等效于pascal代码？

Question

我正在尝试进行一些软件考古。试图用python中的现代脚本替换旧的pascal文件。

在某些时候，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;

在this book, on page 285.中进行了描述，这是范·维恩加登-德克-布伦特方法。

我想只用python中的一行替换它，最好使用scipy。我看到有scipy.optimize.brentq method，但有一个巨大的区别：

• Pascal的zbrent仅使用一个公差输入（tol），而python具有rtol和xtol。我不明白他们的意思。

该怎么办？我应该在我的python程序中给xtol和rtol些什么，使其等同于pascal？我对数值计算一无所知。我很害怕。我只是材料科学家！

Answer 1

Pascal代码中的公差tol作为绝对公差唯一地应用于剩余包围间隔的长度。其功能对应于scipy中的xtol参数。您可以忽略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`.