%% Solve Nonlinear System Without and Including Jacobian
% This example shows the reduction in function evaluations when you provide derivatives for a system of nonlinear equations. As explained in Writing Vector and Matrix Objective Functions, the Jacobian J(x) of a system of equations F(x) is \[{J_{ij}}(x) = \frac{{\partial {F_i}(x)}}{{\partial {x_j}}}\]. Provide this derivative as the second output of your objective function.
% For example, the multirosenbrock function is an n-dimensional generalization of Rosenbrock's function (see Solve a Constrained Nonlinear Problem) for any positive even value of n:
% \[\begin{array}{l}
% F(1) = 1 - {x_1}\\
% F(2) = 10\left( {{x_2} - x_1^2} \right)\\
% F(3) = 1 - {x_3}\\
% F(4) = 10\left( {{x_4} - x_3^2} \right)\\
% \vdots \\
% F(n - 1) = 1 - {x_{n - 1}}\\
% F(n) = 10\left( {{x_n} - x_{n - 1}^2} \right).
% \end{array}\]
% The solution of the equation system F(x) = 0 is the point x_i = 1, i=1\ldots n.
% For this objective function, all of the Jacobian terms J_{ij}(x) are zero except for the terms where i and j differ by at most one. For odd values of i < n, the nonzero terms are
% J_{ii}(x) = -1\\
% J_{(i+1)i}=-20x_i\\
% J_{(i+1)(i+1)}=10.
% The multirosenbrock helper function at the end of this example creates the objective function F(x) and its Jacobian J(x).
% Solve the system of equations starting from the point x_i = -1.9 for odd values of i<n, and x_i = 2 for even values of i. Choose n = 64.

n = 64;  
x0(1:n,1) = -1.9; 
x0(2:2:n,1) = 2;
[x,F,exitflag,output,JAC] = fsolve(@multirosenbrock,x0);

% Examine the distance of the computed solution x from the true solution, and the number of function evaluations that fsolve takes to compute the solution.

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

% fsolve finds the solution, and takes over 1000 function evaluations to do so.
% Solve the system of equations again, this time using the Jacobian. To do so, set the 'SpecifyObjectiveGradient' option to true.

opts = optimoptions('fsolve','SpecifyObjectiveGradient',true);
[x2,F2,exitflag2,output2,JAC2] = fsolve(@multirosenbrock,x0,opts);

% Again, examine the distance of the computed solution x2 from the true solution, and the number of function evaluations that fsolve takes to compute the solution.

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

% fsolve returns the same solution as before, but takes about 20 function evaluations to do so, rather than about 1000. In general, using the Jacobian can save computation and can provide increased robustness, although this example does not show the robustness improvement.
%
% This code creates the multirosenbrock helper function.

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

% Copyright 2020 The MathWorks, Inc.