## Nonlinear Equations with Jacobian

Consider the problem of finding a solution to a system of nonlinear equations whose Jacobian is sparse. The dimension of the problem in this example is 1000. The goal is to find x such that F(x) = 0. Assuming n = 1000, the nonlinear equations are

`$\begin{array}{c}F\left(1\right)=3{x}_{1}-2{x}_{1}^{2}-2{x}_{2}+1,\\ F\left(i\right)=3{x}_{i}-2{x}_{i}^{2}-{x}_{i-1}-2{x}_{i+1}+1,\\ F\left(n\right)=3{x}_{n}-2{x}_{n}^{2}-{x}_{n-1}+1.\end{array}$`

To solve a large nonlinear system of equations, F(x) = 0, you can use the trust-region reflective algorithm available in `fsolve`, a large-scale algorithm (Large-Scale vs. Medium-Scale Algorithms).

### Step 1: Write a file nlsf1.m that computes the objective function values and the Jacobian.

```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```

### Step 2: Call the solve routine for the system of equations.

```xstart = -ones(1000,1); fun = @nlsf1; options = optimoptions(@fsolve,'Display','iter',... 'Algorithm','trust-region',... 'SpecifyObjectiveGradient',true,'PrecondBandWidth',0); [x,fval,exitflag,output] = fsolve(fun,xstart,options);```

A starting point is given as well as the function name. The default method for `fsolve` is trust-region-dogleg, so it is necessary to specify `'Algorithm'` as `'trust-region'` in the `options` argument in order to select the trust-region algorithm. Setting the `Display` option to `'iter'` causes `fsolve` to display the output at each iteration. Setting `'SpecifyObjectiveGradient'` to `true`, causes `fsolve` to use the Jacobian information available in `nlsf1.m`.

The commands display this output:

``` Norm of First-order Iteration Func-count f(x) step optimality 0 1 1011 19 1 2 16.1942 7.91898 2.35 2 3 0.0228027 1.33142 0.291 3 4 0.000103359 0.0433329 0.0201 4 5 7.3792e-07 0.0022606 0.000946 5 6 4.02299e-10 0.000268381 4.12e-05 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.```

A linear system is (approximately) solved in each major iteration using the preconditioned conjugate gradient method. Setting `PrecondBandWidth` to 0 in `options` means a diagonal preconditioner is used. (`PrecondBandWidth` specifies the bandwidth of the preconditioning matrix. A bandwidth of 0 means there is only one diagonal in the matrix.)

From the first-order optimality values, fast linear convergence occurs. The number of conjugate gradient (CG) iterations required per major iteration is low, at most five for a problem of 1000 dimensions, implying that the linear systems are not very difficult to solve in this case (though more work is required as convergence progresses).

If you want to use a tridiagonal preconditioner, i.e., a preconditioning matrix with three diagonals (or bandwidth of one), set `PrecondBandWidth` to the value `1`:

```options = optimoptions(@fsolve,'Display','iter','SpecifyObjectiveGradient',true,... 'Algorithm','trust-region','PrecondBandWidth',1); [x,fval,exitflag,output] = fsolve(fun,xstart,options);```

In this case the output is

``` Norm of First-order Iteration Func-count f(x) step optimality 0 1 1011 19 1 2 16.0839 7.92496 1.92 2 3 0.0458181 1.3279 0.579 3 4 0.000101184 0.0631898 0.0203 4 5 3.16615e-07 0.00273698 0.00079 5 6 9.72481e-10 0.00018111 5.82e-05 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.```

Note that although the same number of iterations takes place, the number of PCG iterations has dropped, so less work is being done per iteration. See Preconditioned Conjugate Gradient Method.

Setting `PrecondBandWidth` to `Inf` (this is the default) means that the solver uses Cholesky factorization rather than PCG.