最近在做中心线提取相关工作,用到求函数零点,看numberical recipes in c上的二分法求零点,思路很简单
原书伪码:
#include <math.h>
#define JMAX 40 Maximum allowed number of bisections.
float rtbis(float (*func)(float), float x1, float x2, float xacc)
Using bisection, nd the root of a function func known to lie between x1 and x2. The root,
returned as rtbis, will be rened until its accuracy is xacc.
{
void nrerror(char error_text[]);
int j;
float dx,f,fmid,xmid,rtb;
f=(*func)(x1);
fmid=(*func)(x2);
if (f*fmid >= 0.0) nrerror("Root must be bracketed for bisection in rtbis");
rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2); Orient the search so that f>0
for (j=1;j<=JMAX;j++) { lies at x+dx.
fmid=(*func)(xmid=rtb+(dx *= 0.5)); Bisection loop.
if (fmid <= 0.0) rtb=xmid;
if (fabs(dx) < xacc || fmid == 0.0) return rtb;
}
nrerror("Too many bisections in rtbis");
return 0.0; Never get here.
}
写了个matlab版本
function result=bitsection(func,x1,x2,xacc) f=func(x1); fmid=func(x2); if(f<0) dx=x2-x1; rtb=x1; else dx=x1-x2; rtb=x2; end for j=1:40 dx=dx*1/2; xmid=rtb+dx; fmid=func(xmid); if(fmid<0) rtb=xmind; end if(abs(dx)<xacc||fmid==0) result=rtb; end end end