2

我有一个黑盒函数 f(x) 和 x 的一系列值。我需要找到 f(x) = 0 的 x 的最小值。

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

我知道 f(x) 是连续的，并且对于所讨论的范围应该只有 0,1 或 2 个根，但它可能有一个局部最小值。

f(x) 在计算上有点昂贵，我必须经常找到这个第一个根。

我在想某种具有一定程度随机性的爬山，以避免任何局部最小值，但是你怎么知道最小值是否小于零，或者你只是还没有找到它？我认为该功能不应该有超过两个最低点，但我不能绝对确定足以依赖它。

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

4

4 回答 4

4

我的方法听起来很复杂，但最终该方法的计算时间将远远小于距离计算（评估你的f(x)）。此外，现有库中已经编写了很多它的实现。

所以我会做什么：

  • f(x)用切比雪夫多项式近似
  • 找到那个多项式的真正根
  • 如果找到任何根，则将这些根用作更精确的根查找器中的初始估计（如果需要）

鉴于您的函数的性质（平滑、连续、否则表现良好）以及存在 0,1 或 2 个根的信息，已经可以通过 3 次评估找到一个好的切比雪夫多项式f(x)。

然后求切比雪夫系数伴生矩阵的特征值；这些对应于切比雪夫多项式的根。

如果一切都是虚构的，则有 0 个根。如果有一些真正的根源，请检查两个是否相等（您所说的“罕见”情况）。否则，所有实特征值都是根；其中最低的一个是你寻找的根。

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

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

我在 Matlab/Octave 中实现了 Chebychev 例程（如下所示）。像这样使用：

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

在[x_min,x_max]您的范围内x，5用于查找多项式的点数（越高，越准确。等于所需的函数评估量），最后一个true将绘制实际函数和切比雪夫近似值的图（主要用于测试目的）。

现在，实现：

% 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
于 2013-05-16T12:07:43.953 回答
0

您可以使用割线法，它是牛顿法的离散版本。

在此处输入图像描述

通过计算两点之间的线（= 割线）及其与 X 轴的交叉点来估计根。

于 2013-05-11T20:16:39.920 回答
0

您的函数只有 0、1 或 2 个根，因此可以使用不能确保第一个根的算法来完成。

  • 使用牛顿法或其他方法求一个根。如果它找不到任何根，这个算法也放弃了。
  • 让找到的根是r并且 x 的范围的开始是x0。让d = (r-x0)/2.
  • 同时d > 0，计算f(r-d)。if f(r-d) > 0, half d( d := d / 2) 和循环。如果 f(r-d) <= 0，退出循环。
  • 如果循环由 完成d = 0，则报告r为第一个根。如果，则使用任何其他方法在和d > 0之间找到一个根并报告它。x0r-d

我假设了两个先决条件。

  1. f(x) 取 x 个浮点数
  2. 在 f(x) 的根的每个点，f(x) 的图形都穿过x 轴。他们没有像 f(x)=x^2 中的 x=0 那样接触根。

使用条件 2，您可以证明如果不存在f(r-d) < 0, ∀ x: x0 < x < r, f(x) > 0。

于 2013-05-12T02:02:03.917 回答
0

您可以uniroot.all对 R 库 rootSolve 中的函数进行小幅更改。

uniroot.all <- function (f, interval, lower= min(interval),
                         upper= max(interval), tol= .Machine$double.eps^0.2,
                         maxiter= 1000, n = 100, nroots = -1, ... ) {

  ## error checking as in uniroot...
  if (!missing(interval) && length(interval) != 2)
    stop("'interval' must be a vector of length 2")
  if (!is.numeric(lower) || !is.numeric(upper) || lower >=
      upper)
    stop("lower < upper  is not fulfilled")

  ## subdivide interval in n subintervals and estimate the function values
  xseq <- seq(lower,upper,len=n+1)
  mod  <- f(xseq,...)

  ## some function values may already be 0
  Equi <- xseq[which(mod==0)]

  ss   <- mod[1:n]*mod[2:(n+1)]  # interval where functionvalues change sign
  ii   <- which(ss<0)

  for (i in ii) {
    Equi <- c(Equi, uniroot(f, lower = xseq[i], upper = xseq[i+1] ,...)$root)
    if (length(Equi) == nroots) {
      return(Equi)
    }
  }
  return(Equi)
}

并像这样运行它：

uniroot.all(f = your_function, interval = c(start, stop), nroots = 1)
于 2016-09-26T17:53:11.373 回答