函数求根次数最少的单变量寻根算法

2024-03-28 21:14:26 发布

您现在位置:Python中文网/ 问答频道 /正文

我正在寻找根查找算法,使用很少的函数评估(目标是最小的)。寻根问题具有以下特征:

f(x) = 0, R -> R

  • 函数(f(.))的评估是非常昂贵的*
  • 边界间隔([a,b])可用于开始(相对较好的近似值,而不是胡乱猜测)
  • f(.)是连续的
  • f(.)是可微的(解析导数不可用)
  • 已知只有一个根位于起始间隔([a,b])内
  • 平滑变化f(.)(函数中不存在极端情况)
  • 允许的停止标准,例如|f(x)| < 1e-2就足够了。你知道吗

*我们可以安全地假设,与f(.)的单个计算相比,该算法所做的任何计算都可以忽略不计。因此,节省甚至一个单一的功能评估是一个重大的收获。你知道吗

考虑到这些,用最少的函数求值找到根的最有效算法是什么?

基于Matlab的^{}和scipy的^{},Brent的方法似乎是流行的选择,尽管对于上面描述的特定问题可能有更有效的算法。你知道吗

也欢迎参考书籍和评论文章。你知道吗


Tags: 函数功能算法目标标准间隔情况scipy
1条回答
网友
1楼 · 发布于 2024-03-28 21:14:26

好吧,如果你想试试布伦特的方法,你可以试试下面我翻译的原始算法。这只是我把原来的C代码翻译成MATLAB。我已经验证了所有原始C代码的测试用例在我的翻译中产生相同的结果。你知道吗

在下面的代码中,This是函数句柄,ab是搜索边界。你知道吗

function x = brent(This,a,b) 

tol = 1e-2 ;
maxit = 1e+3 ;
displayIter = true ;

Fa = This(a) ;
Fb = This(b) ;
c = a ;
Fc = Fa ;

it = 0 ;
if displayIter
    fprintf(1,'Searching for a root in the interval [%g,%g] using Brent''s method...\n',a,b) ;
end
while it<maxit
    it = it + 1 ;

    prevStep = b - a ;

    if abs(Fc) < abs(Fb)
        % swap for b to be best approximation
        a = b ;
        b = c ;
        c = a ;
        Fa = Fb ;
        Fb = Fc ;
        Fc = Fa ;
    end

    tolAct =  2*eps*abs(b) + tol/2 ;

    newStep = (c-b)/2 ;

    if abs(newStep) <= tolAct || abs(Fb)<eps
        % acceptable solution found
        x = b ;
        return 
    end

    if abs(prevStep) >= tolAct && Fa == Fb
        % previous step was large enough and in the right direction, try
        % interpolation
        cb = c-b ;
        if abs(a-c)<eps
            % if only two distinct points, interpolate linearly
            t1 = Fb/Fa ;
            p = cb*t1 ;
            q = 1 - t1 ;
        else
            % three points, do quadratic inverse  interpolation
            a = Fa/Fc ;
            t1 = Fb/Fc ;
            t2 = Fb/Fa ;
            p = t2*( cb*q*(q-t1) - (b-a)*(t1-1) ) ;
            q = (q-1)*(t1-1)*(t2-1) ;
        end
        if p>0
            q = -q ;
        else
            p = -p ;
        end
        if p < ( 0.75*cb*q-abs(tolAct*q)/2 ) && p < abs(prevStep*q/2)
            newStep = p/q ;
        end
    end
    % step must be at least as large as tolerance
    if abs(newStep)<tolAct
        if newStep>0
            newStep = tolAct ;
        else
            newStep = -tolAct ;
        end
    end

    a = b ;
    Fa = Fb ;
    b = b + newStep ;
    Fb = This(b) ;
    if ( Fb > 0 && Fc > 0 ) || ( Fb < 0 && Fc < 0 )
        c = a ;
        Fc = Fa ;
    end
    if displayIter
        fprintf(1,'[%g] Gap = %g, dx = %g\n',it,abs(Fb),newStep) ;
    end
end
fprintf(1,'[%g] Gap = %g, dx = %g\n',it,abs(Fb),newStep) ;

end

相关问题 更多 >