Symbolic Math Toolbox

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.

$$f(\bf{x}) = \sum\limits_{i=1}^{n/2}{100(x_{2i}-x_{2i-1}^2)^2 + (1-x_{2i-1})^2}$$

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

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

  1. Translate the equations for the Rosenbrock function and it's Jacobian from the math form in the text to numeric code.
  2. 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]