【3.18】数值算法-bitsection-二分法求f(x)=0函数零点

 

最近在做中心线提取相关工作，用到求函数零点，看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 rened 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