具有雅可比矩阵的大型稀疏非线性方程组

此示例展示如何使用 fsolve 求解器的功能有效地求解大型稀疏方程组。该示例使用针对 $n$ 方程组定义的目标函数，

$\begin{array}{l} F (1) = 3 x_{1} - 2 x_{1}^{2} - 2 x_{2} + 1, \\ F (i) = 3 x_{i} - 2 x_{i}^{2} - x_{i - 1} - 2 x_{i + 1} + 1, \\ F (n) = 3 x_{n} - 2 x_{n}^{2} - x_{n - 1} + 1 . \end{array}$

要解的方程是 $F_{i} (x) = 0$ 、 $1 \leq i \leq n$ 。该示例使用 $n = 1000$ 。

这个目标函数足够简单，您可以通过分析计算它的雅可比矩阵。如编写向量和矩阵目标函数中所述，方程组 $F (x)$ 的雅可比矩阵 $J (x)$ 是 $J_{i j} (x) = \frac{\partial F_{i} (x)}{\partial x_{j}}$ 。请提供此导数作为目标函数的第二个输出。此示例末尾的 nlsf1 辅助函数会创建目标函数 $F (x)$ 及其雅可比矩阵 $J (x)$ 。

使用默认选项求解方程组，即调用 'trust-region-dogleg' 算法。从点 xstart(i) = -1 开始。

n = 1000;
xstart = -ones(n,1);
fun = @nlsf1; 
[x,fval,exitflag,output] = fsolve(fun,xstart);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

显示解的质量和进行的函数计算次数。

disp(norm(fval))

   2.7603e-13

disp(output.funcCount)

fsolve 可以准确地求解方程，但需要进行数千次函数计算。

使用默认算法和 'trust-region' 算法中的雅可比矩阵求解方程。

options = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,fval2,exitflag2,output2] = fsolve(fun,xstart,options);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

options.Algorithm = 'trust-region';
[x3,fval3,exitflag3,output3] = fsolve(fun,xstart,options);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

disp([norm(fval2),norm(fval3)])

   1.0e-08 *

    0.0000    0.1065

disp([output2.funcCount,output3.funcCount])

     7     5

与不使用雅可比矩阵的情况相比，这两种算法都只占用函数计算次数的一小部分。默认算法比 'trust-region' 算法需要更多的函数计算，但默认算法可以得出更准确的答案。

查看将 'PrecondBandWidth' 选项设置为 1 是否会改变 'trust-region' 答案或效率。此设置导致 fsolve 使用三对角预条件子，这对于该三对角方程组应该有效。

options.PrecondBandWidth = 1;
[x4,fval4,exitflag4,output4] = fsolve(fun,xstart,options);

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

<stopping criteria details>

disp(norm(fval4))

   3.1185e-05

disp(output4.funcCount)

disp(output4.cgiterations)

'PrecondBandWidth' 选项设置导致 fsolve 给出稍微不太准确的答案，以残差范数来衡量。函数计算的数量略有增加，从 5 个增加到 6 个。求解器的解过程包含少于十次共轭梯度迭代。

看看 fsolve 使用对角预条件子后的表现如何。

options.PrecondBandWidth = 0;
[x5,fval5,exitflag5,output5] = fsolve(fun,xstart,options);

Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.

<stopping criteria details>

disp(norm(fval5))

   2.0057e-05

disp(output5.funcCount)

disp(output5.cgiterations)

这次残差范数略低，函数计算的次数保持不变。共轭梯度迭代次数从 8 次增加到 19 次，表明此 'PrecondBandWidth' 设置导致求解器做更多的工作。

使用 'levenberg-marquardt' 算法求解方程。

options = optimoptions('fsolve','SpecifyObjectiveGradient',true,'Algorithm','levenberg-marquardt');
[x6,fval6,exitflag6,output6] = fsolve(fun,xstart,options);

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

disp(norm(fval6))

   5.2545e-15

disp(output6.funcCount)

该算法给出了最低的残差范数，并且使用的函数计算仅比最低次数多一次。

汇总结果。

t = table([norm(fval);norm(fval2);norm(fval3);norm(fval4);norm(fval5);norm(fval6)],...
    [output.funcCount;output2.funcCount;output3.funcCount;output4.funcCount;output5.funcCount;output6.funcCount],...
    'VariableNames',["Residual" "Fevals"],...
    'RowNames',["Default" "Default+Jacobian" "Trust-Region+Jacobian" "Trust-Region+Jacobian,BW=1" "Trust-Region+Jacobian,BW=0" "Levenberg-Marquardt+Jacobian"])

t=6×2 table
                                     Residual     Fevals
                                    __________    ______

    Default                         2.7603e-13     7007 
    Default+Jacobian                2.5891e-13        7 
    Trust-Region+Jacobian           1.0646e-09        5 
    Trust-Region+Jacobian,BW=1      3.1185e-05        6 
    Trust-Region+Jacobian,BW=0      2.0057e-05        6 
    Levenberg-Marquardt+Jacobian    5.2545e-15        6

使用以下代码创建 nlsf1 函数。

function [F,J] = nlsf1(x)
% Evaluate the vector function
n = length(x);
F = zeros(n,1);
i = 2:(n-1);
F(i) = (3-2*x(i)).*x(i)-x(i-1)-2*x(i+1) + 1;
F(n) = (3-2*x(n)).*x(n)-x(n-1) + 1;
F(1) = (3-2*x(1)).*x(1)-2*x(2) + 1;
% Evaluate the Jacobian if nargout > 1
if nargout > 1
   d = -4*x + 3*ones(n,1); D = sparse(1:n,1:n,d,n,n);
   c = -2*ones(n-1,1); C = sparse(1:n-1,2:n,c,n,n);
   e = -ones(n-1,1); E = sparse(2:n,1:n-1,e,n,n);
   J = C + D + E;
end
end

另请参阅

fsolve

具有雅可比矩阵的大型稀疏非线性方程组

另请参阅

主题