# Using Symbolic Math in Numerical Solution of Nonlinear Equations

This example demonstrates that Symbolic Math Toolbox helps minimize errors when solving a nonlinear system of equations. The goal is to find the minimum of the Rosenbrock function, or the so-called "banana" function.

This example uses the following Symbolic Math Toolbox capablities:

- Computing the analytical gradient using
`gradient` - Computing the analytical jacobian using
`jacobian` - Converting symbolic results into numeric functions using
`matlabFunction`

## Contents

- Motivation
- Rewrite Rosenbrock Function as a System of Nonlinear Equations
- Calculate Jacobian of Matrix Representing the System of Equations
- Convert System of Equations and Jacobian to a MATLAB function
- Use fsolve to Numerically Solve the System of Nonlinear Equations
- Use vpasolve to Numerically Solve the System of Nonlinear Equations

## Motivation

The “Nonlinear Equations with Analytic Jacobian” example in Optimization Toolbox solves the same problem by using the `fsolve` function. The workflow shown in that example has two potential sources of errors. You would need to

- Translate the equations for the Rosenbrock function and it's Jacobian from the math form in the text to numeric code.
- Explicitly calculate the Jacobian. For complicated equations, this task is prone to errors.

This example shows that using symbolic math to describe the problem statement and subsequent steps minimizes or even eliminates such errors.

## Rewrite Rosenbrock Function as a System of Nonlinear Equations

First, convert the Rosenbrock function `f` to a system of nonlinear equations `F` , where `F` is a gradient of `f` .

```
n = 64;
x = sym('x',[n,1]);
f = sum(100*(x(2:2:end)-x(1:2:end).^2).^2 + (1-x(1:2:end)).^2);
F = gradient(f,x);
```

Show the first 10 equations:

F(1:10)

ans = 2*x1 - 400*x1*(- x1^2 + x2) - 2 - 200*x1^2 + 200*x2 2*x3 - 400*x3*(- x3^2 + x4) - 2 - 200*x3^2 + 200*x4 2*x5 - 400*x5*(- x5^2 + x6) - 2 - 200*x5^2 + 200*x6 2*x7 - 400*x7*(- x7^2 + x8) - 2 - 200*x7^2 + 200*x8 2*x9 - 400*x9*(- x9^2 + x10) - 2 - 200*x9^2 + 200*x10

To match the intermediate results shown in the Optimization Toolbox example, use the same form of `F`.

F(1:2:end) = simplify(F(1:2:end) + 2*x(1:2:end).*F(2:2:end)); F(1:2:end) = -F(1:2:end)/2; F(2:2:end) = F(2:2:end)/20;

Now the systems of equations are identical:

F(1:10)

ans = 1 - x1 - 10*x1^2 + 10*x2 1 - x3 - 10*x3^2 + 10*x4 1 - x5 - 10*x5^2 + 10*x6 1 - x7 - 10*x7^2 + 10*x8 1 - x9 - 10*x9^2 + 10*x10

## Calculate Jacobian of Matrix Representing the System of Equations

Use `jacobian` to calculate the Jacobian of `F` . This function calculates the Jacobian symbolically, thus avoiding errors associated with numerical approximations of derivatives.

JF = jacobian(F,x);

Show the first 10 rows and columns of the Jacobian matrix.

JF(1:10,1:10)

ans = [ -1, 0, 0, 0, 0, 0, 0, 0, 0, 0] [ -20*x1, 10, 0, 0, 0, 0, 0, 0, 0, 0] [ 0, 0, -1, 0, 0, 0, 0, 0, 0, 0] [ 0, 0, -20*x3, 10, 0, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, -1, 0, 0, 0, 0, 0] [ 0, 0, 0, 0, -20*x5, 10, 0, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, -1, 0, 0, 0] [ 0, 0, 0, 0, 0, 0, -20*x7, 10, 0, 0] [ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0] [ 0, 0, 0, 0, 0, 0, 0, 0, -20*x9, 10]

Most of the elements of the Jacobian matrix `JF` are zeros. Nevertheless, when you convert this matrix to a MATLAB function, the result is a dense numeric matrix. Often, operations on sparse matrices are more efficient than the same operations on dense matrices.

Therefore, to create efficient code, convert only the nonzero elements of `JF` to a symbolic vector before invoking matlabFunction. `is` and `js` represent the sparsity pattern of `JF`.

[is,js] = find(JF); JF = JF(JF~=0);

## Convert System of Equations and Jacobian to a MATLAB function

The system of equations `F` , representing the Rosenbrock function, is a symbolic matrix that consists of symbolic expressions. To be able to solve it with the `fsolve` function, convert this system to a MATLAB function. At this step, it is convenient to convert both `F` and its Jacobian, `JF` , to a single file-based MATLAB function, `FJFfun`. In principle, this allows `F` and `JF` to reuse variables. Alternatively, you can convert them to individual MATLAB functions.

matlabFunction([F;JF],'var',x,'file','FJFfun');

`FJFfun` accepts variables as a list, but `fsolve` expects a vector of variables. The function `genericObj` converts the vector to a list. Anonymous function `bananaObj` is defined to fix argument values for `genericObj`. Note the use of the sparsity pattern in the function `genericObj`. Further note that this approach of using `FJFfun` and `genericObj` together is not tied to the particular expression represented by `F` and can be used as is for any `F`.

```
bananaobj = @(y) genericObj(y,@FJFfun,is,js);
type genericObj
```

% genericObj takes as input, vector x, Function and Jacobian evaluation % function handle FJFfun, and sparsity pattern is,js and returns as output, % the numeric values of the Function and Jacobian: F and J function [F,J] = genericObj(x,FJFfun,is,js) xcell = num2cell(x); % FJs(1:length(x)) is the numeric vector for the Function % FJs(length(x)+1:end) is the numeric vector for the non-zero entries of the Jacobian FJs = FJFfun(xcell{:}); F = FJs(1:length(x)); J = sparse(is,js,FJs(length(x)+1:end)); end

## Use fsolve to Numerically Solve the System of Nonlinear Equations

Use `fsolve` for the system of equations, converted to MATLAB function. Use the starting point *x(i) = –1.9* for the odd indices, and *x(i) = 2* for the even indices. Set `'Display'` to `'iter'` to see the solver's progress. Set `'Jacobian'` to `'on'` to use the Jacobian defined in `bananaobj`.

x0(1:n,1) = -1.9; x0(2:2:n,1) = 2; options = optimoptions(@fsolve,'Display','iter','Jacobian','on'); [sol1,Fsol,exitflag,output,JAC] = fsolve(bananaobj,x0,options);

Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 1 8563.84 615 1 1 2 3093.71 1 329 1 2 3 225.104 2.5 34.8 2.5 3 4 212.48 6.25 34.1 6.25 4 5 212.48 6.25 34.1 6.25 5 6 102.771 1.5625 6.39 1.56 6 7 102.771 3.90625 6.39 3.91 7 8 87.7443 0.976563 2.19 0.977 8 9 74.1426 2.44141 6.27 2.44 9 10 74.1426 2.44141 6.27 2.44 10 11 52.497 0.610352 1.52 0.61 11 12 41.3297 1.52588 4.63 1.53 12 13 34.5115 1.52588 6.97 1.53 13 14 16.9716 1.52588 4.69 1.53 14 15 8.16797 1.52588 3.77 1.53 15 16 3.55178 1.52588 3.56 1.53 16 17 1.38476 1.52588 3.31 1.53 17 18 0.219553 1.16206 1.66 1.53 18 19 0 0.0468565 0 1.53 Equation solved. fsolve completed because the vector of function values is near zero as measured by the default value of the function tolerance, and the problem appears regular as measured by the gradient.

## Use vpasolve to Numerically Solve the System of Nonlinear Equations

The system of nonlinear equations, `F`, can be alternatively solved with the `vpasolve` function, which is a numeric solver available in Symbolic Math Toolbox. A key difference between `vpasolve` and `fsolve` is that `fsolve` uses double-precision arithmetic to solve the equations, whereas `vpasolve` uses variable-precision arithmetic and therefore, lets you control the accuracy of computations.

sol2 = vpasolve(F);

If you assign the solution of a system of equations to one variable, `sol2`, then `vpasolve` returns the solutions in form of a structure array. You can access each solution (each field of the structure array):

sol2.x1

ans = 1.0

You also can convert `sol2` to a symbolic vector, and then access a range of the solutions. For example, display the 5th to 25th solutions:

for k=(1:64) sol2_array(k) = sol2.(char(x(k))); end sol2_array(5:25)

ans = [ 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]