确定fsolve的收敛速度

Question

假设我正在解一个非线性方程组。一个简单的例子是：

function example
    x0 = [15; -2];
    options = optimoptions('fsolve','Display','iter','TolFun',eps,'TolX',eps);
    [x,fval,exitflag,output] = fsolve(@P1a,x0,options);
end

function f1 = P1a(x)
    f1 = [x(1)+x(2)*(x(2)*(5-x(2))-2)- 13; x(1)+x(2)*(x(2)*(1+x(2))-14)-29];
end

我如何确定收敛速度呢？'Display'，'iter'在每一步都向我展示了规范，但我找不到提取这些值的方法。(对于这个特殊的例子，我认为fsolve没有收敛到正确的解决方案，而是收敛到一个局部最小值。然而，这不是问题所在。我只想找到一种估计收敛速度的方法。)

Answer 1

你可以从fsolve中学到很多东西。但是，您需要做一些工作。阅读'OutputFcn'选项并为Matlab的优化方法编写输出函数。这与Matlab的ODE求解器使用的相同名称的选项相似。下面是一个示例，它复制了'Display','iter'选项- fsolve的值(特别是用于默认的'trust-region-dogleg'算法)：

function stop = outfun(x,optimValues,state)
% See private/trustnleqn
stop = false;

switch state
    case 'init'
        header = sprintf(['\n                                         Norm of      First-order   Trust-region\n',...
                          ' Iteration  Func-count     f(x)          step         optimality    radius']);
        disp(header);
    case 'iter'
        iter = optimValues.iteration;               % Iteration
        numFevals = optimValues.funccount;          % Func-count
        F = optimValues.fval;                       % f(x)
        normd = optimValues.stepsize;               % Norm of step
        normgradinf = optimValues.firstorderopt;    % First-order optimality
        Delta = optimValues.trustregionradius;      % Trust-region radius

        if iter > 0
            formatstr = ' %5.0f      %5.0f   %13.6g  %13.6g   %12.3g    %12.3g';
            iterOutput = sprintf(formatstr,iter,numFevals,F'*F,normd,normgradinf,Delta);
        else
            formatstr0 = ' %5.0f      %5.0f   %13.6g                  %12.3g    %12.3g';
            iterOutput = sprintf(formatstr0,iter,numFevals,F'*F,normgradinf,Delta);
        end
        disp(iterOutput);
    case 'done'

    otherwise

end

然后，您可以通过以下方式进行调用：

function example
P1a=@(x)[x(1)+x(2)*(x(2)*(5-x(2))-2)- 13; x(1)+x(2)*(x(2)*(1+x(2))-14)-29];
x0 = [15; -2];
opts = optimoptions('fsolve','Display','off','OutputFcn',@outfun,'TolFun',eps,'TolX',eps);
[x,fval,exitflag,output] = fsolve(P1a,x0,opts);

这仍然只是打印到命令窗口。从这里开始，创建一个可以将数据写入数组、文件或其他数据结构的输出函数。下面是如何使用全局变量 (通常情况下，不是一个好主意)：

function stop = outfun2(x,optimValues,state)
stop = false;
global out;   % Global variable, define in main function too

switch state
    case 'init'
        out = [];
    case 'iter'
        iter = optimValues.iteration;               % Iteration
        numFevals = optimValues.funccount;          % Func-count
        F = optimValues.fval;                       % f(x)
        normd = optimValues.stepsize;               % Norm of step
        normgradinf = optimValues.firstorderopt;    % First-order optimality
        Delta = optimValues.trustregionradius;      % Trust-region radius

        out = [out;iter numFevals F'*F normd normgradinf Delta];
    case 'done'

    otherwise

end

然后在调用global out;之前先在主函数中声明fsolve。您还可以通过将输出函数设为嵌套函数来实现这一点，在这种情况下，out数组将与外部主函数共享。

第二个输出函数示例执行动态内存分配，而不是重新分配整个out数组。这是无法避免的，因为我们和算法都不知道要进行多少迭代才能收敛。然而，在几百次迭代中，动态内存分配将是非常快的。

我把“决定收敛速度”留给你们，既然你们已经掌握了这些工具.