时间戳

首页 > 解决方案 > 如何使用 scipy.brentq 函数，使其等效于帕斯卡代码？

python - 如何使用 scipy.brentq 函数，使其等效于帕斯卡代码？

问题描述

我正在尝试做一些软件考古。试图用 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;

本书第 285 页对此进行了描述。它是 Van Wijngaarden-Dekker-Brent 方法。

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

该怎么办？我应该在我的 python 程序中给出什么作为 xtol 和 rtol 以使其等效于帕斯卡？我对数值计算一无所知。我很害怕。我只是一个材料科学家！

标签： pythonscipynumerical-methods

解决方案

帕斯卡代码中的公差tol唯一地适用于作为绝对公差的剩余括号间隔的长度。它的功能对应xtolscipy中的参数。您可以忽略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`.

推荐阅读