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))
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;
14. fc=fb;
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;
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. }
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

明白了，谢谢您。