# Preface

Nonlinear System is a system in which the change of the output is not proportional to the change of the input, which can modeled with a set of nonlinear equations. Solving the nonlinear equations can give us the clues of the behavior of a nonlinear system. Most of the time, the system is so complex that we can not solve it analytically but only numerically.

Assume there is a simple system of nonlinear equation:

And we need to solve it numerically. $\alpha$, $\beta$, $\lambda$, $\gamma$ are parameters.

The system of nonlinear equations can be written into myfun.m as follow:

To input parameters $\alpha$, $\beta$, $\lambda$, $\gamma$ into myfun, we can use para vector or global command in MATLAB. $\alpha$, $\beta$, $\lambda$, $\gamma$ are numerical data. When there are string type of data in parameters, we can use cell instead of vector.

If $\alpha \neq 0$ and $\lambda \neq 0$, then we can get four solutions as follow: $(x,y)=(-1,5)$, $(x,y)=(-1,-7)$, $(x,y)=(3,5)$, $(x,y)=(3,-7)$. Which one to be chosen depends on the initial condition. Thus initial condition is really important for solving system of nonlinear equations.

In most cases, the system of nonlinear equations that we would like to solve is really complicated that we can not even guess how the analytical solution could be. For example, if we want to solve $\epsilon^d$, $\epsilon^c$, $\theta$ and $alpha$ from the system of nonlinear equations as follow:

Parameters and functions $F(x)$ and $q (\theta)$ are stated as follow:

The functions $u$, $h_s$ and $h_n$ are defined as:

It should be emphasized that we also regard $\sigma$ as a variable. ($\sigma \in [0,0.5]$) Our final goal, we want to figure out the dynamic relationship of $\epsilon^d$, $\epsilon^c$, $\theta$, $alpha$ $\sim$ $\sigma$.

We are going to show how to solve this problem via using MATLAB, Fortran, Python and R.

# MATLAB: fsolve, fmincon and lsqnonlin

In MATLAB, there are 3 kinds of solver to solve nonlinear equations: fsolve, fmincon and lsqnonlin.0

var vector represent $[\epsilon^d;\epsilon^c;\theta;\alpha]$.
parameters vector represents $[c^F;c^p;\beta;\phi;\delta;\sigma;\lambda;b;r;\epsilon^c;\mu;step;width;typen]$.
We can define the nonlinear equation set as as follow:

And we can also define function $F(x)$, $q(\theta)$, $u$, $h_s$ and $h_n$ as follow:

## fsolve

Solves a problem specified by F(x) = 0 for x, where F(x) is a function that returns a vector value. x is a vector or a matrix; see Matrix Arguments.

fsolve is the most commonly used solver for nonlinear equations in MATLAB. In MATLAB, we can find three kinds of algorithms: Trust Region Dogleg1, Trust region and Levenberg-Marquardt.

We can set the options of solver fsolve as follow:

The nonlinear equations that we want to solve are coded in function handle F.

And then we can solve nonlinear equation:

where fs is the vector of solution and ceq is the value of nonlinear functions.

## fmincon

Find minimum of constrained nonlinear multivariable function

For fmincon, there are 3 kinds of algorithms: Interior Point2, Trust Region Reflective3, Sequential Quadratic Programming (SQP)45, Active Set.

We can set the options as follow:

We also need to set some constraint vectors following the suggested format of fmincon.

And then we can solve nonlinear equation:

where fun_solveJun01T1_nonlcon.m is as follow:

## lsqnonlin

Solve nonlinear least-squares (nonlinear data-fitting) problems

In lsqnonlin solver, Trust Region Reflective3 and Levenberg-Marquardt are used.

We can set the options as follow:

options = optimoptions(@lsqnonlin);
options.Algorithm = 'trust-region-reflective';
str_algo = 'TRR';
%   'trust-region-reflective' (default)
%   A gradient to be supplied in the objective function
%   'levenberg-marquardt'
options.MaxFunctionEvaluations = Inf;
options.MaxIterations = 500;

We can also use function handle F to represent the nonlinear equations:

F = @(var)fun_solveJun01T1(var,para);
% var = [epsilon_d;epsilon_c;theta]
lblsq = [-Inf;-Inf;-Inf;-Inf];
if strcmp(typen,'I')
ublsq = [epsilon_u;epsilon_u;Inf;Inf];
else
ublsq = [epsilon_u;epsilon_u;Inf;Inf];
end

where lblsq and ublsp are the lower and upper bound.

## Results

We can solve the equation Jun01T1 and get the solution for $\epsilon^d$, $\epsilon^c$, $\theta$ and $alpha$. In main.m:

Then we can get plot the solution as follow

i. For the fsolve solver:

ii. For the fmincon solver:

iii. For the lsqnonlin solver:

# Fortran: NEQNF

In IMSL Numerical Library, NEQNF can be used to solve a system of nonlinear equations by Powell’s Method, which uses finite-difference approximation to approximate Jacobian Matrix.

## NEQNF

We can define the functions $u$, $h_s$ and $h_n$ in subroutines fun_u, fun_hs and fun_hn as follow:

In subroutine Jun01T1, we can use CALL DNEQNF(FCN,ERRREL,N,ITMAX,XGUESS,X,fnorm) to solve system of nonlinear equations FCN.

where subroutine FCN is defined as:

In main.f90, we call the subroutine Jun01T1, fun_u, fun_hs and fun_hn:

## Results

Then we can plot the solutions as follow:

# Python: fsolve in scipy.optimize package

We can use fsolve in scipy.optimize, which is a wrapper around [MINPACK’s hybrd and hybrj algorithms.

The purpose of HYBRD is to find a zero of a system of N non-linear functions in N variables by a modification of the Powell hybrid method. The user must provide a subroutine which calculates the functions. The Jacobian is then calculated by a forward-difference approximation.

The purpose of HYBRJ is to find a zero of a system of N non-linear functions in N variables by a modification of the Powell hybrid method. The user must provide a subroutine which calculates the functions and the Jacobian.

HYBRD is a modification of the Powell hybrid method. Two of its main characteristics involve the choice of the correction as a convex combination of the Newton and scaled gradient directions, and the updating of the Jacobian by the rank-1 method of Broyden.

## fsolve in scipy.optimize

We define a Class named Operator to solve the nonlinear equations:

In the function Solution(self), we use fsolve in the format as follow:

## Results

We can take use of matplotlib.pyplot to plot the solutions as follow:

# R: nleqslv package

To solve system of nonlinear equations, we can use nleqslv package. The nleqslv package provides two algorithms for solving (dense) nonlinear systems of equations:

• a Broyden Secant method6 where the matrix of derivatives is updated after each major iteration using the Broyden rank 1 update.
• a full Newton method where the Jacobian matrix of derivatives is recalculated at each iteration

## nleqslv package

We can define the nonlinear Equations to be solved as follow:

To solve the nonlinear equations, we can use nleqslv as follow:

## Results

Taking use of the plot function in R, we can plot the solutions as follow:

# Mathematica: Export results as image or .m file

Generally, we can use NSolve in Mathematica to solve system of nonlinear equations. However, it really takes time to handle with the nonlinear equations that we want to solve. But there are a few useful techniques that I would like to cover here.

## Export result as images

We can define the function $F(x)$ and $q(\theta)$ as follow:

We can also define the system of nonlinear equations and variables as follow:

Then we can calculate the Jacobian Matrix as follow:

We can export the results of partial derivatives as images:

## Export as MATLAB file

Or we can take use of the Wolfram Library ToMatlab to convert the results into MATLAB .m files as follow:

For example, dF_depsilond.m file is as follow

ToMatlab in Mathematica can really make MATLAB much more powerful because of the better analytical capacity of Mathematica compared with MATLAB.

## Results

The exported image of partial derivatives, such as $\frac{\partial F}{\partial \epsilon^d}$, is as follow:

We can plot the dynamics between partial derivatives and some variables (such as $\alpha$) as follow:

# More: To ODE and PDE

System of Nonlinear Equation has big relationship with Ordinary Differential Equations (ODE) and Partial Differential Equations (PDE) in numerical analysis. And there are so many famous nonlinear equations out there in various disciplines.7 Thus knowing how to solve nonlinear equations numerically is crucial in mathematical modeling in different aspects.

0. MathWork:Nonlinear Systems with Constraints
1. Lecture 7: The Dogleg and Steihaug Methods
2. Interior-point method for NLP
3. Trust Region Reflective Algorithm
4. Constrained Nonlinear Optimization Algorithms
5. Nocedal, J., & Wright, S. (2006). Numerical optimization. Springer Science & Business Media.
6. Broyden’s method
7. Some famous nonlinear equations
-------------End of postThanks for your time-------------
Enjoy it? Subscribe to my blog by scanning my public wechat account
Rate this post