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

shingyuan

我最近在研究指数函数方程的数值解法，我看的数是"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 )}}$
这个值作为新的对根的近似估计，如果出现一些不好的情况，则使用二分法代替以上插值法

但该算法实现的时候，出现了这样的情况

1. s=fb/fa;
   
2.                         if(a==c)
   
3.                         {
   
4.                                 p=2.0*xm*s;
   
5.                                 q=1.0-s;
   
6.                         }
   
7.                         else
   
8.                         {
   
9.                                 q=fa/fc;
   
10.                                 r=fb/fc;
    
11.                                 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
    
12.                                 q=(q-1.0)*(r-1.0)*(s-1.0);
    
13.                         }

复制代码

其中

1. xm=0.5*(c-b);

复制代码

书上所给的代码，计算 P/Q的时候，P差了一个负号，相比较于公式。这是什么原因？

这里的判断也不是很好理解

1. if(p>0.0 )
   
2. {
   
3.          q= -q;
   
4.   }

复制代码

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


下面有两句不太好理解的代码：

1. min1=3.0*xm*q-fabs(tol1*q);
   
2.                         min2=fabs(e*q);

复制代码

这两个极值算出来是什么意义？

下面接着的这个条件判断

1. 2.0*p<(min1 < min2 ? min1 : min2)

复制代码

似乎也不太好理解，请教下大家，谢谢

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

附：算法的源代码

1. #include <math.h>
   
2. #include <stdio.h>
   
3. #define ITMAX        100
   
4. #define EPS 3.0e-8
   
5. #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
   
6. 
   
7. double zbrent(double (*func)(double), double x1, double x2, double tol)
   
8. {
   
9.         int iter;
   
10.         double a=x1,b=x2,c=x2,d,e,min1,min2;
    
11.         double fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm;
    
12. 
    
13.        
    
14.         fc=fb;
    
15.        
    
16.         for (iter=1;iter<=ITMAX;iter++)
    
17.         {
    
18.                 if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0))
    
19.                 {
    
20.                                 c=a;
    
21.                                 fc=fa;
    
22.                                 e=d=b-a;
    
23.                 }
    
24.                 if (fabs(fc) < fabs(fb))
    
25.                 {
    
26.                                 a=b;
    
27.                                 b=c;
    
28.                                 c=a;
    
29.                                 fa=fb;
    
30.                                 fb=fc;
    
31.                                 fc=fa;
    
32.                 }
    
33.                 tol1=2.0*EPS*fabs(b)+0.5*tol;
    
34.                 xm=0.5*(c-b);
    
35.                 if(fabs(xm)<=tol1 || fb==0.0) return b;
    
36.                
    
37.                 if(fabs(e) >= tol1 && fabs(fa) > fabs(fb))
    
38.                 {
    
39.                         s=fb/fa;
    
40.                         if(a==c)
    
41.                         {
    
42.                                 p=2.0*xm*s;
    
43.                                 q=1.0-s;
    
44.                         }
    
45.                         else
    
46.                         {
    
47.                                 q=fa/fc;
    
48.                                 r=fb/fc;
    
49.                                 p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
    
50.                                 q=(q-1.0)*(r-1.0)*(s-1.0);
    
51.                         }
    
52.                         if(p>0.0)
    
53.                         {
    
54.                                 q= -q;
    
55.                         }
    
56. 
    
57.                         p=fabs(p);
    
58.                         min1=3.0*xm*q-fabs(tol1*q);
    
59.                         min2=fabs(e*q);
    
60.                         if(2.0*p<(min1 < min2 ? min1 : min2))
    
61.                         {
    
62.                                 e=d;
    
63.                                 d=p/q;
    
64.                         }
    
65.                         else
    
66.                         {
    
67.                                 d=xm;
    
68.                                 e=d;
    
69.                         }
    
70.                 }
    
71.                 else
    
72.                 {
    
73.                         d=xm;
    
74.                         e=d;
    
75.                 }
    
76.                 a=b;
    
77.                 fa=fb;
    
78.                 if(fabs(d)>tol1)
    
79.                 {
    
80.                         b+=d;
    
81.                 }
    
82.                 else
    
83.                 {
    
84.                         b+=SIGN(tol1,xm);
    
85.                 }
    
86.                 fb=(*func)(b);
    
87.         }
    
88.         return 0.0;
    
89. }

复制代码

gxqcn

...
附：为何 \left( \right)  \left[ \right] 这些会没有用？？
shingyuan 发表于 2011-4-26 13:51

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

shingyuan

Web 上的 LaTeX 是简易版，为的是保证低带宽下满足基本的交流。
gxqcn 发表于 2011-4-26 14:59

明白了，谢谢您。