查找黑盒函数的第一根，或同一函数的任何负值 [英] Find first root of a black box function, or any negative value of same function

问题描述

我有一个黑盒函数f(x)和x的取值范围. 我需要找到f(x)= 0的x的最小值.

I have a black box function, f(x) and a range of values for x. I need to find the lowest value of x for which f(x) = 0.

我知道对于x范围的开始，f(x)> 0，并且如果我有一个f(x)< 0我可以使用regula falsi或类似的寻根方法来尝试确定f(x)= 0.

I know that for the start of the range of x, f(x) > 0, and if I had a value for which f(x) < 0 I could use regula falsi, or similar root finding methods, to try determine f(x)=0.

我知道f(x)是连续的，并且该范围的根仅应有0,1或2个根，但是它可能具有局部最小值.

I know f(x) is continuous, and should only have 0,1 or 2 roots for the range in question, but it might have a local minimum.

f(x)在计算上有些昂贵，我将不得不找到很多第一个根.

f(x) is somewhat computationally expensive, and I'll have to find this first root a lot.

我当时正在考虑某种程度的随机性爬山，以避免任何局部最小值，但是您如何知道最小值是否小于零，或者您是否尚未找到最小值?我认为该函数不应有两个以上的最低点，但是我不能绝对肯定要依赖它.

I was thinking some kind of hill climbing with a degree of randomness to avoid any local minimums, but then how do you know if there was no minimum less than zero or if you just haven't found it yet? I think the function shouldn't have more than two minimum points, but I can't be absolutely certain of that enough to rely on it.

如果有帮助，在这种情况下，x表示时间，而f(x)表示该时间船与轨道上的物体(月球/行星)之间的距离.我需要第一个点是它们彼此之间有一定距离.

If it helps, x in this case represents a time, and f(x) represents the distance between a ship and a body in orbit (moon/planet) at that time. I need the first point where they are a certain distance from each other.

推荐答案

我的方法听起来很复杂，但最终该方法的计算时间将比距离计算小得多(对f(x)的求值) .此外，现有库中已经编写了很多实现.

My method will sound pretty complicated, but in the end the computation time of the method will be far smaller than the distance calculations (evaluation of your f(x)). Also, there are quite many implementations of it already written up in existing libraries.

那我该怎么做:

• 具有切比雪夫多项式的近似f(x)
• 找到该多项式的实根
• 如果找到任何根，请使用这些根作为更精确的根查找器中的初始估计值(如果需要)

鉴于函数的性质(平滑，连续，否则行为良好)以及存在0、1或2个根的信息，可以通过3个对f(x)的求值来找到一个好的Chebychev多项式.

Given the nature of your function (smooth, continuous, otherwise well-behaved) and the information that there's 0,1 or 2 roots, a good Chebychev polynomial can already be found with 3 evaluations of f(x).

然后找到切比雪夫系数的伴随矩阵的特征值；这些对应于切比雪夫多项式的根.

Then find the eigenvalues of the companion matrix of the Chebychev coefficients; these correspond to the roots of the Chebychev polynomial.

如果全部都是虚构的，则存在0个根. 如果存在真正的根，请检查两个是否相等(您所说的稀有"情况). 否则，所有真实特征值都是根.最低的是您要寻找的根.

If all are imaginary, there's 0 roots. If there are some real roots, check if two are equal (that "rare" case you spoke of). Otherwise, all real eigenvalues are roots; the lowest one of which is the root you seek.

然后使用Newton-Raphson进行细化(如有必要，或使用更好的Chebychev多项式). f的导数可以使用中心差来近似

Then use Newton-Raphson to refine (if necessary, or use a better Chebychev polynomial). Derivatives of f can be approximated using central differences

f'(x) = ( f(x+h)-f(h-x) ) /2/h     (for small h)

我在Matlab/Octave中实现了Chebychev例程(如下所示).像这样使用:

I have an implementation of the Chebychev routines in Matlab/Octave (given below). Use like this:

R = FindRealRoots(@f, x_min, x_max, 5, true,true);

，其中[x_min,x_max]是您在x中的范围，5用于查找多项式的点数(越高，越精确.等于所需的函数求值数量)，最后一个true将绘制实际功能及其切比雪夫近似值(主要用于测试).

with [x_min,x_max] your range in x, 5 the number of points to use for finding the polynomial (the higher, the more accurate. Equals the amount of function evaluations needed), and the last true will make a plot of the actual function and the Chebychev approximation to it (mainly for testing purposes).

现在，实现:

% FINDREALROOTS     Find approximations to all real roots of any function 
%                   on an interval [a, b].
%
% USAGE:
%    Roots = FindRealRoots(funfcn, a, b, n, vectorized, make_plot)
%
% FINDREALROOTS() approximates all the real roots of the function 'funfcn' 
% in the interval [a,b]. It does so by finding the roots of an [n]-th degree 
% Chebyshev polynomial approximation, via the eignevalues of the associated 
% companion matrix. 
%
% When the argument [vectorized] is [true], FINDREALROOTS() will evaluate 
% the function 'funfcn' at all [n] required points in the interval 
% simultaneously. Otherwise, it will use ARRAFUN() to calculate the [n] 
% function values one-by-one. [vectorized] defaults to [false]. 
%
% When the argument [make_plot] is true, FINDREALROOTS() plots the 
% original function and the Chebyshev approximation, and shows any roots on
% the given interval. Also [make_plot] defaults to [false]. 
%
% All [Roots] (if any) will be sorted.
% 
% First version 26th May 2007 by Stephen Morris, 
% Nightingale-EOS Ltd., St. Asaph, Wales.
%
% Modified 14/Nov (Rody Oldenhuis)
%
% See also roots, eig.

function Roots = FindRealRoots(funfcn, a, b, n, vectorized, make_plot)

    % parse input and initialize.
    inarg = nargin; 
    if n <= 2, n = 3; end                    % Minimum [n] is 3:    
    if (inarg < 5), vectorized = false; end  % default: function isn't vectorized
    if (inarg < 6), make_plot = false; end   % default: don't make plot

    % some convenient variables
    bma = (b-a)/2;  bpa = (b+a)/2;  Roots = [];

    % Obtain the Chebyshev coefficients for the function
    %
    % Based on the routine given in Numerical Recipes (3rd) section 5.8;
    % calculates the Chebyshev coefficients necessary to approximate some
    % function over the interval [a,b]

    % initialize 
    c = zeros(1,n);  k=(1:n)';  y = cos(pi*((1:n)-1/2)/n); 
    % evaluate function on Chebychev nodes
    if vectorized
        f = feval(funfcn,(y*bma)+bpa);
    else
        f = arrayfun(@(x) feval(funfcn,x),(y*bma)+bpa);
    end

    % compute the coefficients
    for j=1:n, c(j)=(f(:).'*(cos((pi*(j-1))*((k-0.5)/n))))*(2-(j==1))/n; end       

    % coefficients may be [NaN] if [inf]
    % ??? TODO - it is of course possible for c(n) to be zero...
    if any(~isfinite(c(:))) || (c(n) == 0), return; end

    % Define [A] as the Frobenius-Chebyshev companion matrix. This is based
    % on the form given by J.P. Boyd, Appl. Num. Math. 56 pp.1077-1091 (2006).
    one = ones(n-3,1);
    A = diag([one/2; 0],-1) + diag([1; one/2],+1);
    A(end, :) = -c(1:n-1)/2/c(n);
    A(end,end-1) = A(end,end-1) + 0.5;

    % Now we have the companion matrix, we can find its eigenvalues using the
    % MATLAB built-in function. We're only interested in the real elements of
    % the matrix:
    eigvals = eig(A);  realvals = eigvals(imag(eigvals)==0);

    % if there aren't any real roots, return
    if isempty(realvals), return; end

    % Of course these are the roots scaled to the canonical interval [-1,1]. We
    % need to map them back onto the interval [a, b]; we widen the interval just
    % a tiny bit to make sure that we don't miss any that are right on the 
    % boundaries.
    rangevals = nonzeros(realvals(abs(realvals) <= 1+1e-5));

    % also sort the roots
    Roots = sort(rangevals*bma + bpa);

    % As a sanity check we'll plot out the original function and its Chebyshev
    % approximation: if they don't match then we know to call the routine again
    % with a larger 'n'.
    if make_plot
        % simple grid
        grid = linspace(a,b, max(25,n));
        % evaluate function
        if vectorized
            fungrid = feval(funfcn, grid);
        else
            fungrid = arrayfun(@(x) feval(funfcn,x), grid);
        end        
        % corresponding Chebychev-grid (more complicated but essentially the same)
        y = (2.*grid-a-b)./(b-a); d = zeros(1,length(grid)); dd = d;
        for j = length(c):-1:2, sv=d; d=(2*y.*d)-dd+c(j); dd=sv; end, chebgrid=(y.*d)-dd+c(1);
        % Now make plot
        figure(1), clf,  hold on
        plot(grid, fungrid ,'color' , 'r');
        line(grid, chebgrid,'color' , 'b'); 
        line(grid, zeros(1,length(grid)), 'linestyle','--')
        legend('function', 'interpolation')
    end % make plot

end % FindRealRoots

这篇关于查找黑盒函数的第一根，或同一函数的任何负值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！