Main Content

求解不含和含雅可比矩阵的非线性方程组

此示例说明在为非线性方程组提供导数时函数计算次数会减少。如编写向量和矩阵目标函数中所述,方程组 F(x) 的雅可比矩阵 J(x)Jij(x)=Fi(x)xj。请提供此导数作为目标函数的第二个输出。

例如,multirosenbrock 函数是罗森布罗克函数的 n 维泛化(请参阅基于问题求解有约束非线性问题:),适用于 n 的任何正偶数值:

F(1)=1-x1F(2)=10(x2-x12)F(3)=1-x3F(4)=10(x4-x32)F(n-1)=1-xn-1F(n)=10(xn-xn-12).

方程组 F(x)=0 的解是点 xi=1i=1n

对于此目标函数,除了 ij 最多相差一之外,所有雅可比项 Jij(x) 都为零。对于 i<n 的奇数值,非零项为

Jii(x)=-1J(i+1)i=-20xiJ(i+1)(i+1)=10.

此示例末尾multirosenbrock 辅助函数会创建目标函数 F(x) 及其雅可比矩阵 J(x)

对于 i<n 的奇数值,从点 xi=-1.9 开始求解方程组;对于 i 的偶数值,从点 xi=2 开始求解方程组。指定 n=64

n = 64;  
x0(1:n,1) = -1.9; 
x0(2:2:n,1) = 2;
[x,F,exitflag,output,JAC] = fsolve(@multirosenbrock,x0);
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.

检查计算的解 x 与实际解的距离,以及 fsolve 计算该解所执行的函数计算次数。

disp(norm(x-ones(size(x))))
     0
disp(output.funcCount)
        1043

fsolve 找到了解,并为此执行了 1000 多次函数计算。

再次求解方程组,这次使用雅可比矩阵。为此,请将 'SpecifyObjectiveGradient' 选项设置为 true

opts = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,F2,exitflag2,output2,JAC2] = fsolve(@multirosenbrock,x0,opts);
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.

同样,检查计算的解 x2 与实际解的距离,以及 fsolve 计算该解所执行的函数计算次数。

disp(norm(x2-ones(size(x2))))
     0
disp(output2.funcCount)
    21

fsolve 返回与上一个解相同的解,但为此执行的函数计算是 20 次左右,而不是 1000 多次。一般情况下,使用雅可比矩阵可以减少函数计算次数并提高稳健性,尽管此示例并未说明稳健性得到了改进。

辅助函数

以下代码会创建 multirosenbrock 辅助函数。

function [F,J] = multirosenbrock(x)
% Get the problem size
n = length(x);  
if n == 0, error('Input vector, x, is empty.'); end
if mod(n,2) ~= 0
   error('Input vector, x ,must have an even number of components.');
end
% Evaluate the vector function
odds  = 1:2:n;
evens = 2:2:n;
F = zeros(n,1);
F(odds,1)  = 1-x(odds);
F(evens,1) = 10.*(x(evens)-x(odds).^2); 
% Evaluate the Jacobian matrix if nargout > 1
if nargout > 1
   c = -ones(n/2,1);    C = sparse(odds,odds,c,n,n);
   d = 10*ones(n/2,1);  D = sparse(evens,evens,d,n,n);
   e = -20.*x(odds);    E = sparse(evens,odds,e,n,n);
   J = C + D + E;
end
end

另请参阅

相关主题