论坛 › 算法交流 › Brent算法的一些细节问题(公式编辑好了)

shingyuan 发表于 2011-4-26 13:51:57

Brent算法的一些细节问题(公式编辑好了)

本帖最后由 shingyuan 于 2011-4-26 15:30 编辑

我最近在研究指数函数方程的数值解法，我看的数是"Numerical Recipes in C", 第九章介绍了相关方法，我试验的结果是"Bisection"方法速度稍慢，"Section"方法有的时候区间会超出原始区间，带来一些问题。
Brent方法好像没有什么太大问题，所以我决定使用这种方法，这个算法的想法应该是对三个初始点进行二次插值，其零点解出为：
x=b+\frac{\frac{f(b)}{f(a)}(\frac{f(a)}{f(c)}(\frac{f(b)}{f(c)}-\frac{f(a)}{f(c)})(c-b)-(1-\frac{f(b)}{f(c)})(b-a))}{(\frac{f(a)}{f(c)}-1)(\frac{f(b)}{f(c)}-1)(\frac{f(b)}{f(a)}-1)}


对于有两个点相同的情形为：x=b+\frac{\frac{f ( b )}{f ( a )}( b-a)}{1-\frac{f ( b )}{f ( a t )}}
这个值作为新的对根的近似估计，如果出现一些不好的情况，则使用二分法代替以上插值法

但该算法实现的时候，出现了这样的情况s=fb/fa;
                        if(a==c)
                        {
                                p=2.0*xm*s;
                                q=1.0-s;
                        }
                        else
                        {
                                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);
                        }其中xm=0.5*(c-b);书上所给的代码，计算 P/Q的时候，P差了一个负号，相比较于公式。这是什么原因？

这里的判断也不是很好理解 if(p>0.0 )
{
         q= -q;
}这里按照书上的说法是b经过小的修正P/Q后（P是x = b + */*的分子,Q是分母）,是否会越过初始的区间,这样做是否是判断P/Q的符号？另外，判断后为啥将q变号？是否是补上前面缺少的一个负号，感觉不太好理解


下面有两句不太好理解的代码：min1=3.0*xm*q-fabs(tol1*q);
                        min2=fabs(e*q);这两个极值算出来是什么意义？

下面接着的这个条件判断2.0*p<(min1 < min2 ? min1 : min2)似乎也不太好理解，请教下大家，谢谢

附：为何 \left( \right)\left[ \right] 这些会没有用？？

附：算法的源代码#include <math.h>
#include <stdio.h>
#define ITMAX        100
#define EPS 3.0e-8
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

double zbrent(double (*func)(double), double x1, double x2, double tol)
{
        int iter;
        double a=x1,b=x2,c=x2,d,e,min1,min2;
        double fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm;

       
        fc=fb;
       
        for (iter=1;iter<=ITMAX;iter++)
        {
                if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0))
                {
                                c=a;
                                fc=fa;
                                e=d=b-a;
                }
                if (fabs(fc) < fabs(fb))
                {
                                a=b;
                                b=c;
                                c=a;
                                fa=fb;
                                fb=fc;
                                fc=fa;
                }
                tol1=2.0*EPS*fabs(b)+0.5*tol;
                xm=0.5*(c-b);
                if(fabs(xm)<=tol1 || fb==0.0) return b;
               
                if(fabs(e) >= tol1 && fabs(fa) > fabs(fb))
                {
                        s=fb/fa;
                        if(a==c)
                        {
                                p=2.0*xm*s;
                                q=1.0-s;
                        }
                        else
                        {
                                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);
                        }
                        if(p>0.0)
                        {
                                q= -q;
                        }

                        p=fabs(p);
                        min1=3.0*xm*q-fabs(tol1*q);
                        min2=fabs(e*q);
                        if(2.0*p<(min1 < min2 ? min1 : min2))
                        {
                                e=d;
                                d=p/q;
                        }
                        else
                        {
                                d=xm;
                                e=d;
                        }
                }
                else
                {
                        d=xm;
                        e=d;
                }
                a=b;
                fa=fb;
                if(fabs(d)>tol1)
                {
                        b+=d;
                }
                else
                {
                        b+=SIGN(tol1,xm);
                }
                fb=(*func)(b);
        }
        return 0.0;
}

gxqcn 发表于 2011-4-26 14:59:06

...
附：为何 \left( \right)\left[ \right] 这些会没有用？？
shingyuan 发表于 2011-4-26 13:51 http://bbs.emath.ac.cn/images/common/back.gif

Web 上的 LaTeX 是简易版，为的是保证低带宽下满足基本的交流。

shingyuan 发表于 2011-4-26 15:01:30

Web 上的 LaTeX 是简易版，为的是保证低带宽下满足基本的交流。
gxqcn 发表于 2011-4-26 14:59 http://bbs.emath.ac.cn/images/common/back.gif

明白了，谢谢您。

查看完整版本: Brent算法的一些细节问题(公式编辑好了)