Solving System of Nonlinear Equations: using MATLAB, Fortran, Python and R language

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:

myfun.m
1
2
3
4
5
6
7
8
9
10
function F = myfun(var,para)
alpha=para(1);
beta=para(2);
lambda=para(3);
gamma=para(4);
x=var(1);
y=var(2);
F(1) = alpha*(x+1)*(x-3)*exp(beta*y);
F(2) = lambda*(y-5)*(y+7)*exp(gamma*x);
return

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.

Click to see Source Code

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:

fun_solveJun01T1.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
function F = fun_solveJun01T1(var,parameters)
dbstop if error
A = parameters{1};
B1= parameters{2};
B2 = parameters{3};
c_f = parameters{4};
c_p = parameters{5};
beta = parameters{6};
phi = parameters{7};
delta = parameters{8};
sigma = parameters{9};
lambda = parameters{10};
b = parameters{11};
r = parameters{12};
epsilon_u = parameters{13};
mu = parameters{14};
pstar = parameters{15};
typen = parameters{16};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon_d=var(1);
epsilon_c=var(2);
theta=var(3);
alpha=var(4);
p=pstar;
q_theta = fun_q_theta(theta,A,B1,B2);
int_Fdu = fun_int_F(epsilon_d,epsilon_u,typen,epsilon_u);
int_Fcu = fun_int_F(epsilon_c,epsilon_u,typen,epsilon_u);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
part11 = epsilon_d+lambda*int_Fdu/(r+lambda+theta*q_theta);
part12 = (b-p)/sigma;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
temp211 = delta*(r+lambda-phi*(1-fun_F_x(epsilon_c,typen,epsilon_u)))...
*(epsilon_c-epsilon_d);
temp212 = r+lambda+theta*q_theta;
temp213 = (lambda/(r+lambda)-lambda*delta/temp212)*int_Fcu;
part21 = epsilon_c-temp211/((1-phi)*temp212)+temp213;
temp2211 = ((beta+phi*(1-beta))*c_f*(1-alpha))/((1-beta)*(1-phi));
temp2212 = (beta*c_p*alpha)/(1-beta);
temp221 = theta*(temp2211 + temp2212);
part22 = delta*epsilon_d+((1-delta)*(b-p)+temp221)/delta;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
part31 = 1/q_theta;
temp321 = (1-beta)*(1-phi)/c_f;
temp322 = sigma*(epsilon_u-epsilon_c)/(r+lambda)+...
delta*sigma*(epsilon_c-epsilon_d)/((1-phi)*(r+lambda+theta*q_theta));
part32 = temp321*temp322;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
part41 = 1/q_theta;
temp421 = (1-beta)*delta*sigma*(epsilon_u-epsilon_d);
temp422 = c_p*(r+lambda+theta*q_theta);
part42 = temp421/temp422;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = [part11-part12;part21-part22;part31 - part32;part41-part42];
return

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

fun_F_x.m
1
F_x(i)=.5-.5*erf((log(-x(i)+1)+.5*log(2))/sqrt(2*log(2)));
fun_q_theta.m
1
2
3
4
5
function q_theta = fun_q_theta(theta,A,B1,B2)
temp1 = nthroot(theta,B2);
temp2 = (temp1).^(-B1);
q_theta = A*temp2;
return
fun_u.m
1
2
3
4
5
6
function u = fun_u(epsilon_d,theta,lambda,typen,...
A,B1,B2,epsilon_u)
q_theta=fun_q_theta(theta,A,B1,B2);
F_epsilond = fun_F_x(epsilon_d,typen,epsilon_u);
u = lambda*F_epsilond/(lambda*F_epsilond+theta*q_theta);
return
fun_hs.m
1
2
3
4
5
6
7
8
9
function h_s = fun_hs(epsilon_c,epsilon_d,theta,alpha,u,lambda,typen,...
A,B1,B2,epsilon_u)
q_theta=fun_q_theta(theta,A,B1,B2);
F_epsilonc = fun_F_x(epsilon_c,typen,epsilon_u);
F_epsilond = fun_F_x(epsilon_d,typen,epsilon_u);
temp1 = (1-u)*(1-alpha)*lambda*(F_epsilonc-F_epsilond);
temp2 = lambda*F_epsilonc+(1-alpha)*theta*q_theta;
h_s =temp1 / temp2;
return
fun_hn.m
1
2
3
4
5
6
7
8
function h_n = fun_hn(epsilon_c,theta,alpha,u,lambda,typen,...
A,B1,B2,epsilon_u)
q_theta=fun_q_theta(theta,A,B1,B2);
F_epsilonc = fun_F_x(epsilon_c,typen,epsilon_u);
temp1 = (1-u)*(alpha)*lambda*(F_epsilonc);
temp2 = lambda*F_epsilonc+(1-alpha)*theta*q_theta;
h_n =temp1 / temp2;
return

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:

Jun01T1.m
1
2
3
4
5
6
7
8
options = optimoptions('fsolve');
options.Algorithm = 'trust-region-dogleg';
% 'trust-region-dogleg' (default)
% 'trust-region-reflective'
% 'levenberg-marquardt'
options.MaxFunEvals = Inf;
options.MaxIter = 500;
options.Display = 'iter';

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

Jun01T1.m
1
2
F = @(var)fun_solveJun01T1(var,para);
% var = [epsilon_d;epsilon_c;theta]

And then we can solve nonlinear equation:

Jun01T1.m
1
2
[fs,~,exitflag] = fsolve(F,var0,options);
ceq=fun_solveJun01T1(fs,para);

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:

Jun01T1.m
1
2
3
4
5
options = optimoptions(@fmincon);
options.Algorithm = 'interior-point';
str_algo = 'IP';
options.MaxFunctionEvaluations = Inf;
options.MaxIterations = 500;

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

Jun01T1.m
1
2
3
4
5
6
7
Afmin=[];bfmin=[];Aeqfmin=[];beqfmin=[];
lbfmin = [-Inf;-Inf;-Inf;-Inf];
if strcmp(typen,'I')
ubfmin = [epsilon_u;epsilon_u;Inf;Inf];
else
ubfmin = [epsilon_u;epsilon_u;Inf;Inf];
end

And then we can solve nonlinear equation:

Jun01T1.m
1
2
3
[fs,~,exitflag] = fmincon(@(var)0,var0,Afmin,bfmin,Aeqfmin,beqfmin,lbfmin,ubfmin,...
@(x)fun_solveJun01T1_nonlcon(x,para),options);
[~,ceq]=fun_solveJun01T1_nonlcon(fs,para);

where fun_solveJun01T1_nonlcon.m is as follow:

fun_solveJun01T1_nonlcon.m
1
2
3
4
5
function [c,ceq] = fun_solveJun01T1_nonlcon(var,para)
c=[]; % no nonlinear inequality
ceq = fun_solveJun01T1(var,para);
% the fsolve objective is fmincon constrain
return

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:

main.m
1
2
3
4
5
6
7
8
9
10
11
12
para ={A;B1;B2;c_f;c_p;beta;phi;delta;sigma;lambda;...
b;r;epsilon_u;mu;pstar;typen;taskcode;...
str_s;str_para0;path_fig;vt;i};
[vepsilon_d(i),vepsilon_c(i),vtheta(i),valpha(i),vceq(:,i),...
vexitflag(i),str_algo]=Jun01T1(para);
str_para=strcat(str_para0,'-',str_algo);
vu(i) = fun_u(vepsilon_d(i),vtheta(i),lambda,typen,...
A,B1,B2,epsilon_u);
vhs(i)=fun_hs(vepsilon_c(i),vepsilon_d(i),vtheta(i),valpha(i),...
vu(i),lambda,typen,A,B1,B2,epsilon_u);
vhn(i)=fun_hn(vepsilon_c(i),vtheta(i),valpha(i),vu(i),...
lambda,typen,A,B1,B2,epsilon_u);

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:

fun_u.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
subroutine fun_u(epsilon_d,theta,u)
USE point_para
IMPLICIT NONE
REAL(kind=8),intent(in) :: epsilon_d,theta
REAL(kind=8),intent(out) :: u
REAL(kind=8),external :: funqtheta,fun_F_x
! Local variables
REAL(kind=8) :: q_theta,F_epsilond,A,B1,B2,lambda
!!!!!!!!!!!!!
A = p_para(1)
B1 = p_para(14)
B2 = p_para(15)
lambda = p_para(9)
q_theta = funqtheta(theta,A,B1,B2)
F_epsilond = fun_F_x(epsilon_d)
u = lambda*F_epsilond/(lambda*F_epsilond+theta*q_theta)
return
end subroutine
fun_hs.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
subroutine fun_hs(epsilon_c,epsilon_d,theta,alpha,u,h_s)
USE point_para
implicit none
REAL(kind=8),intent(in) :: epsilon_c,epsilon_d,theta,alpha,u
REAL(kind=8),intent(out) :: h_s
REAL(kind=8),external :: funqtheta,fun_F_x
! LOCAL VARIABLES
REAL(kind=8) :: q_theta,F_epsilonc,F_epsilond,temp1,temp2,A,B1,B2,lambda
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
A = p_para(1)
B1 = p_para(14)
B2 = p_para(15)
lambda = p_para(9)
!!!!!!!!!!!!!!!!!!!!!!!!!!!
q_theta = funqtheta(theta,A,B1,B2)
F_epsilonc = fun_F_x(epsilon_c)
F_epsilond = fun_F_x(epsilon_d)
temp1 = (1-u)*(1-alpha)*lambda*(F_epsilonc-F_epsilond)
temp2 = lambda*F_epsilonc+(1-alpha)*theta*q_theta
h_s =temp1/temp2
return
end subroutine
fun_hn.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
subroutine fun_hn(epsilon_c,theta,alpha,u,h_n)
USE point_para
IMPLICIT NONE
!!!!!!!!!!!!!!!!!!!!!!!!!!
REAL(kind=8),intent(in) :: epsilon_c,theta,alpha,u
REAL(kind=8),intent(out) :: h_n
REAL(kind=8),external :: funqtheta,fun_F_x
! LOCAL VARIABLES
REAL(kind=8) :: q_theta,F_epsilonc,temp1,temp2,A,B1,B2,lambda
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
A = p_para(1)
B1 = p_para(14)
B2 = p_para(15)
lambda = p_para(9)
!!!!!!!!!!!!!!!!!!!!!!!!!!!
q_theta = funqtheta(theta,A,B1,B2)
F_epsilonc = fun_F_x(epsilon_c)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
temp1 = (1-u)*(alpha)*lambda*(F_epsilonc)
temp2 = lambda*F_epsilonc+(1-alpha)*theta*q_theta
h_n =temp1/temp2
return
end subroutine

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

Jun01T1.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
subroutine Jun01T1(epsilond,epsilonc,theta,alpha,exitflag)
INCLUDE 'link_fnl_shared.h'
USE NEQNF_INT
implicit none
real(kind=8) :: epsilond
real(kind=8) :: epsilonc
real(kind=8) :: theta
real(kind=8) :: alpha
logical(kind=4) :: exitflag
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! The number of equations to be solved. for NEQNF
integer(kind=4),parameter :: N=4
real(kind=8):: fnorm, X(N), XGUESS(N)
! fnorm: A scalar that has the value F(1)2 + … + F(N)2 at the point X. (Output of NEQNF)
! X(N): Solution of nonlinear systems; A vector of length N. (Output of NEQNF)
! X contains the best estimate of the root found by NEQNF.
! XGUESS: A vector of length N. (Input of NEQNF)
real(kind=8),parameter:: ERRREL = 1d-6
! ERRREL:Tolerance; Stopping criterion; (Input of NEQNF)
integer(kind=4),PARAMETER:: ITMAX = 1500
! The maximum allowable number of iterations. (Input)
EXTERNAL FCN
! User-supplied SUBROUTINE to evaluate the system of equations to be solved.
! The usage is CALL fun_solveJun01T1(X, F, N)
!real ceq(N)
! A vector of length N. FVEC contains the functions evaluated at the point X.
DATA XGUESS /-5.0,-1.0,1.0,0.5/
! Routine NEQNF is based on the MINPACK subroutine HYBRD1,
! which uses a modification of M.J.D. Powell’s hybrid algorithm.
CALL DNEQNF(FCN,ERRREL,N,ITMAX,XGUESS,X,fnorm)
write(*,*) X
epsilond = X(1)
epsilonc = X(2)
theta = X(3)
alpha = X(4)
! Check whether the solution is valid
if (fnorm<=ERRREL) then
exitflag = .true.
else
exitflag = .false.
end if
return
end subroutine Jun01T1

where subroutine FCN is defined as:

fun_solveJun01T1.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
subroutine FCN (X,F,N)
use point_para
implicit none
integer(KIND=4) :: N
real(kind=8):: X(N),F(N)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(kind=8) :: A,B1,B2,r,c_p,beta,phi,delta,sigma,lambda,pstar,b,c_f,epsilon_u,mu
real(kind=8) :: p
real(kind=8) :: q_theta,int_Fdu,int_Fcu
real(kind=8) :: part11,part12,temp211,temp212,temp213,part21,temp2211,temp2212,temp221,part22,&
part31,part32,temp321,temp322,part41,part42,temp421,temp422
character(len=3) :: typen
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(kind=8) :: fun_intF,fun_F_X,funqtheta
external :: funqtheta,fun_intF,fun_F_X
REAL EXP, SIN
INTRINSIC EXP, SIN
A = p_para(1)
B1= p_para(14)
B2 = p_para(15)
r = p_para(3)
c_p = p_para(4)
beta = p_para(5)
phi = p_para(6)
delta = p_para(7)
sigma = p_para(8)
lambda = p_para(9)
pstar = p_para(10)
b = p_para(11)
c_f = p_para(12)
mu = p_para(13)
typen = p_typen
epsilon_u=p_epu
p=pstar
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
q_theta = funqtheta(X(3),A,B1,B2)
int_Fdu = fun_intF(X(1),epsilon_u)
int_Fcu = fun_intF(X(2),epsilon_u)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
part11 = X(1)+lambda*int_Fdu/(r+lambda+X(3)*q_theta)
part12 = (b-p)/sigma
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
temp211 = delta*(r+lambda-phi*(1-fun_F_x(X(2))))*(X(2)-X(1))
temp212 = r+lambda+X(3)*q_theta
temp213 = (lambda/(r+lambda)-lambda*delta/temp212)*int_Fcu
part21 = X(2)-temp211/((1-phi)*temp212)+temp213
temp2211 = ((beta+phi*(1-beta))*c_f*(1-X(4)))/((1-beta)*(1-phi))
temp2212 = (beta*c_p*X(4))/(1-beta)
temp221 = X(3)*(temp2211 + temp2212)
part22 = delta*X(1)+((1-delta)*(b-p)+temp221)/delta
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
part31 = 1/q_theta
temp321 = (1-beta)*(1-phi)/c_f
temp322 = sigma*(epsilon_u-X(2))/(r+lambda)+&
delta*sigma*(X(2)-X(1))/((1-phi)*(r+lambda+X(3)*q_theta))
part32 = temp321*temp322
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
part41 = 1/q_theta
temp421 = (1-beta)*delta*sigma*(epsilon_u-X(1))
temp422 = c_p*(r+lambda+X(3)*q_theta)
part42 = temp421/temp422
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
F(1) = part11-part12
F(2) = part21-part22
F(3) = part31-part32
F(4) = part41-part42
return
end SUBROUTINE fcn

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

main.f90
1
2
3
4
call Jun01T1(vepsilon_d(i),vepsilon_c(i),vtheta(i),valpha(i),vexitflag(i))
call fun_u(vepsilon_d(i),vtheta(i),vu(i))
call fun_hs(vepsilon_c(i),vepsilon_d(i),vtheta(i),valpha(i),vu(i),vhs(i))
call fun_hn(vepsilon_c(i),vtheta(i),valpha(i),vu(i),vhn(i))

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:

main_Jun01T1.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
class Operator(CaseCode):
"""Main operator to get the solution
epsilon_d,epsilon_c,theta,alpha"""
# para is the dictionary updated
def __init__(self, case, typen, para):
CaseCode.__init__(self, case)
self.typen = typen
self.para = para
def fun_q_theta(self, x):
return self.para['A']* x **(-self.para['B1']/self.para['B2'])
def fun_F_x(self, x):
if self.typen == 'O':
temp1 = (np.pi/4.0)/np.sin(np.pi/4.0) - x
temp2 = (np.pi/2.0)/np.sin(np.pi/2.0) - ((np.pi/4.0)/np.sin(np.pi/4.0))**2.0
temp3 = 1.0 + (temp1/np.sqrt(temp2))**(-4.0)
return 1.0-(temp3)**(-1.0)
elif self.typen == 'I':
if x < -1.0*self.para['epsilon_u']:
return 0
elif x > self.para['epsilon_u']:
return 1.0
else:
return (x + self.para['epsilon_u'])/(2.0*self.para['epsilon_u'])
elif self.typen == 'II':
temp=np.sqrt(1.0/np.pi)-x*np.sqrt(.5-1.0/np.pi)
return 1.0-special.erf(temp)
elif self.typen == 'III':
return (.5-.5 * special.erf((np.log(-x+1.0)+ \
.5 * np.log(2.0))/np.sqrt(2.0 * np.log(2.0))))
elif self.typen == 'IV':
return ((3.0-2.0*x)/np.sqrt(3.0))**(-3.0)
def integrand(self,x):
return (1-self.fun_F_x(x))
def fun_int_F(self,a,b):
temp = integrate.quad(self.integrand,a,b)
return temp[0]
def Equation(self,x):
epsilon_d = x[0]
epsilon_c = x[1]
theta = x[2]
alpha = x[3]
q_theta = self.fun_q_theta(theta)
int_Fdu = self.fun_int_F(epsilon_d,self.para['epsilon_u'])
int_Fcu = self.fun_int_F(epsilon_c,self.para['epsilon_u'])
####
part11 = epsilon_d + self.para['lambda'] * \
int_Fdu / (self.para['r'] + self.para['lambda'] + theta * q_theta)
part12 = (self.para['b'] - self.para['pstar'])/self.para['sigma']
####
temp211 = self.para['delta'] * (self.para['r'] + self.para['lambda']\
-self.para['phi'] * (1.0-self.fun_F_x(epsilon_c)))\
* (epsilon_c-epsilon_d)
temp212 = self.para['r'] + self.para['lambda'] + theta * q_theta
temp213 = (self.para['lambda']/(self.para['r']+self.para['lambda'])\
-self.para['lambda'] * self.para['delta']/temp212) * int_Fcu
part21 = epsilon_c-temp211/((1.0-self.para['phi']) * temp212)+temp213
temp2211 = ((self.para['beta']+self.para['phi'] * \
(1.0-self.para['beta'])) * self.para['c_f'] * (1.0-alpha))/((1.0-self.para['beta']) * \
(1.0-self.para['phi']))
temp2212 = (self.para['beta'] * self.para['c_p'] * alpha)/(1.0-self.para['beta'])
temp221 = theta * (temp2211 + temp2212)
part22 = self.para['delta'] * epsilon_d+((1.0-self.para['delta']) * \
(self.para['b']-self.para['pstar'])+temp221)/self.para['delta']
####
part31 = 1.0/q_theta
temp321 = (1.0-self.para['beta']) * (1.0-self.para['phi'])/self.para['c_f']
temp322 = self.para['sigma'] * (self.para['epsilon_u']-epsilon_c)/(\
self.para['r']+self.para['lambda'])+\
self.para['delta'] * self.para['sigma'] * (epsilon_c-epsilon_d)/(\
(1.0-self.para['phi']) * (self.para['r']+self.para['lambda']+theta * q_theta))
part32 = temp321 * temp322
####
part41 = 1.0/q_theta
temp421 = (1.0-self.para['beta']) * self.para['delta'] * self.para['sigma'] * (\
self.para['epsilon_u']-epsilon_d)
temp422 = self.para['c_p'] * (self.para['r']+self.para['lambda']+theta * q_theta)
part42 = temp421/temp422
return [
part11-part12,
part21-part22,
part31-part32,
part41-part42
]
def Solution(self):
# Tuple
# Initial value guess
guess = (-5.0,-1.0,1.0,0.5)
# fsolve is a wrapper around MINPACK's hybrd and hybrj algorithms
# [epsilon_d,epsilon_c,theta,alpha]
# xtol : The calculation will terminate if the relative error between
# two consecutive iterates is at most xtol.
# maxfev : The maximum number of calls to the function.
sol, infodict, ier, mesg = fsolve(self.Equation,guess,full_output=True)
if ier != 1:
print("ier = {}".format(ier))
print("sol = {}".format(sol))
print("infodict = {}".format(infodict))
return [sol,ier]
else:
return [sol,ier]
def Tell(self):
CaseCode.Tell(self)
print('Type of F(x) is:{}'.format(self.typen))
print('Parameters are:{}'.format(self.para))
print('Solution is:{}'.format(self.Solution))
print('Error is:{}'.format(self.Error))

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

1
sol, infodict, ier, mesg = fsolve(self.Equation,guess,full_output=True)

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:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
Equation<-function(x){
y<-numeric(length(x))
epsilon_d <- x[1]
epsilon_c <- x[2]
theta <- x[3]
alpha <- x[4]
q_theta = fun_q_theta(myPara=myPara,x=theta)
int_Fdu = fun_int_F(epsilon_d,epsilon_u)
int_Fcu = fun_int_F(epsilon_c,epsilon_u)
part11 = epsilon_d+myPara["lambda","value"]*int_Fdu/(myPara["r","value"]+myPara["lambda","value"]+theta*q_theta)
part12 = (myPara["b","value"]-myPara["pstar","value"])/myPara["sigma","value"]
temp211 = myPara['delta',"value"]*(myPara['r',"value"]+myPara['lambda',"value"]-myPara['phi',"value"]*(1-fun_F_x(epsilon_c)))*(epsilon_c-epsilon_d)
temp212 = myPara['r',"value"]+myPara['lambda',"value"]+theta*q_theta
temp213 = (myPara['lambda',"value"]/(myPara['r',"value"]+myPara['lambda',"value"])-myPara['lambda',"value"]*myPara['delta',"value"]/temp212)*int_Fcu
part21 = epsilon_c-temp211/((1-myPara['phi',"value"])*temp212)+temp213
temp2211 = ((myPara['beta',"value"]+myPara['phi',"value"]*(1.0-myPara['beta',"value"]))*myPara['c_f',"value"]*(1.0-alpha))/((1-myPara['beta',"value"])*(1-myPara['phi',"value"]))
temp2212 = (myPara['beta',"value"]*myPara['c_p',"value"]*alpha)/(1-myPara['beta',"value"])
temp221 = theta*(temp2211 + temp2212)
part22 = myPara['delta',"value"]*epsilon_d+((1-myPara['delta',"value"])*(myPara['b',"value"]-myPara['pstar',"value"])+temp221)/myPara['delta',"value"]
part31 = 1/q_theta
temp321 = (1-myPara['beta',"value"])*(1-myPara['phi',"value"])/myPara['c_f',"value"]
temp322 = myPara['sigma',"value"]*(epsilon_u-epsilon_c)/(myPara['r',"value"]+myPara['lambda',"value"])+myPara['delta',"value"]*myPara['sigma',"value"]*(epsilon_c-epsilon_d)/((1-myPara['phi',"value"])*(myPara['r',"value"]+myPara['lambda',"value"]+theta*q_theta))
part32 = temp321*temp322
part41 = 1/q_theta
temp421 = (1-myPara['beta',"value"])*myPara['delta',"value"]*myPara['sigma',"value"]*(epsilon_u-epsilon_d)
temp422 = myPara['c_p',"value"]*(myPara['r',"value"]+myPara['lambda',"value"]+theta*q_theta)
part42 = temp421/temp422
y[1]<-part11-part12
y[2]<-part21-part22
y[3]<-part31-part32
y[4]<-part41-part42
y
}

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

1
2
3
4
5
6
sol<-nleqslv(xstart,Equation,method="Broyden",control = list(trace=1))
epsilon_d[i]<-sol$x[1]
epsilon_c[i]<-sol$x[2]
theta[i]<-sol$x[3]
alpha[i]<-sol$x[4]

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:

1
2
3
4
5
6
7
8
9
10
11
funqtheta[\[Theta]_, A_, B_] := A*\[Theta]^(-B);
funFx[x_, typen_, epsilonu_] :=
Switch[typen,
"III", .5 - .5*Erf[(Log[-x + 1] + .5*Log[2])/(Sqrt[2*Log[2]])]];
funintF[a_, b_, typen_, epsilonu_] :=
Switch[typen, "III",
Integrate[(-Sqrt[2*Log[2]])*(.5 + .5*Erf[y])*
E^(Sqrt[2*Log[2]]*y - Log[2]/2),
{y, (Log[-a + 1] + .5*Log[2])/(Sqrt[
2*Log[2]]), (Log[-b + 1] + .5*Log[2])/(Sqrt[2*Log[2]])}]];
(* Set y = (Log[-x+1]+Log[2]/2)/Sqrt[2*Log[2]] *)

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
part11 = \[Epsilon]d + \[Lambda]*
intFdu/(r + \[Lambda] + \[Theta]*q\[Theta]);
part12 = (b - pstar)/\[Sigma];
(*******************************************************************************)
temp211 = \[Delta]*(r + \[Lambda] - \[Phi]*(1 -
funFx[\[Epsilon]c,
typen, \[Epsilon]u]))*(\[Epsilon]c - \[Epsilon]d);
temp212 = r + \[Lambda] + \[Theta]*q\[Theta];
temp213 = (\[Lambda]/(r + \[Lambda]) - \[Lambda]*\[Delta]/temp212)*
intFcu;
part21 = \[Epsilon]c - temp211/((1 - \[Phi])*temp212) + temp213;
temp2211 = ((\[Beta] + \[Phi]*(1 - \[Beta]))*
cF*(1 - \[Alpha]))/((1 - \[Beta])*(1 - \[Phi]));
temp2212 = (\[Beta]*cP*\[Alpha])/(1 - \[Beta]);
temp221 = \[Theta]*(temp2211 + temp2212);
part22 = \[Delta]*\[Epsilon]d + ((1 - \[Delta])*(b - pstar) +
temp221)/\[Delta];
(*******************************************************************************)
part31 = 1/q\[Theta];
temp321 = (1 - \[Beta])*(1 - \[Phi])/cF;
temp322 = \[Sigma]*(\[Epsilon]u - \[Epsilon]c)/(r + \[Lambda]) + \
\[Delta]*\[Sigma]*(\[Epsilon]c - \[Epsilon]d)/((1 - \[Phi])*(r + \
\[Lambda] + \[Theta]*q\[Theta]));
part32 = temp321*temp322;
(*******************************************************************************)
part41 = 1/q\[Theta];
temp421 = (1 - \[Beta])*\[Delta]*\[Sigma]*(\[Epsilon]u - \[Epsilon]d);
temp422 = cP*(r + \[Lambda] + \[Theta]*q\[Theta]);
part42 = temp421/temp422;
(*******************************************************************************)
F = {part11 - part12, part21 - part22, part31 - part32,
part41 - part42};
var = {\[Epsilon]d, \[Epsilon]c, \[Theta], \[Alpha],
pstar, b, \[Phi], \[Sigma], \[Beta], \[Lambda], cF, r,
cP, \[Delta]};

Then we can calculate the Jacobian Matrix as follow:

1
2
(* Jacobian Matrix *)
Jac = D[F, {var}]

We can export the results of partial derivatives as images:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Simplify[ArrayReshape[Jac[[All, 1]], {4, 1}]] // TraditionalForm
Export["dF_depsilond.gif",
Simplify[ArrayReshape[Jac[[All, 1]], {4, 1}]] // TraditionalForm]
Export["dF_depsilond.png",
Simplify[ArrayReshape[Jac[[All, 1]], {4, 1}]] // TraditionalForm]
Simplify[ArrayReshape[Jac[[All, 2]], {4, 1}]] // TraditionalForm
Export["dF_depsilonc.gif",
Simplify[ArrayReshape[Jac[[All, 2]], {4, 1}]] // TraditionalForm]
Export["dF_depsilonc.png",
Simplify[ArrayReshape[Jac[[All, 2]], {4, 1}]] // TraditionalForm]
Simplify[ArrayReshape[Jac[[All, 3]], {4, 1}]] // TraditionalForm
Export["dF_dtheta.gif",
Simplify[ArrayReshape[Jac[[All, 3]], {4, 1}]] // TraditionalForm]
Export["dF_dtheta.png",
Simplify[ArrayReshape[Jac[[All, 3]], {4, 1}]] // TraditionalForm]
Simplify[ArrayReshape[Jac[[All, 4]], {4, 1}]] // TraditionalForm
Export["dF_dalpha.gif",
Simplify[ArrayReshape[Jac[[All, 4]], {4, 1}]] // TraditionalForm]
Export["dF_dalpha.png",
Simplify[ArrayReshape[Jac[[All, 4]], {4, 1}]] // TraditionalForm]

Export as MATLAB file

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

1
2
3
4
WriteMatlab[Jac[[All, 1]], "dF_depsilond.m", "dF_depsilond", 50]
WriteMatlab[Jac[[All, 2]], "dF_depsilonc.m", "dF_depsilonc", 50]
WriteMatlab[Jac[[All, 3]], "dF_dtheta.m", "dF_dtheta", 50]
WriteMatlab[Jac[[All, 4]], "dF_dalpha.m", "dF_dalpha", 50]

For example, dF_depsilond.m file is as follow

dF_depsilond_fun.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
function dF_depsilond=dF_depsilond_fun(A,B1,B2,mu,epsilon_u,...
c_f,c_p,beta,phi,delta,sigma,lambda,b,r,pstar,...
epsilon_d,epsilon_c,theta,alpha)
% This function was generated by ToMatlab Package from Mathematica
% 18-Jun-2017
B=B1./B2;
dF_depsilond=[1+lambda.*(r+A.*theta.^(1+(-1).*B)+lambda).^(-1).*((-0.5E0)+( ...
-0.479179E0).*exp(1).^((-1).*(0.294353E0+ ...
0.849322E0.*log(0.1E1+(-0.1E1).*epsilon_d)).^2)+ ...
0.479179E0.*exp(1).^((-1).*(0.294353E0+( ...
-0.849322E0).*log(0.1E1+(-0.1E1).*epsilon_d)).^2).*( ...
0.1E1+(-0.1E1).*epsilon_d).^(-1)+(-0.5E0).*erf( ...
0.294353E0+0.849322E0.*log(0.1E1+(-0.1E1).*epsilon_d))),( ...
-1).*delta+delta.*(r+A.*theta.^(1+(-1).*B)+lambda).^(-1).*(1+(-1).* ...
phi).^(-1).*(r+lambda+(-1).*phi.*(0.5E0+0.5E0.*erf((2.*log( ...
2)).^(-1/2).*(0.346574E0+log(1+(-1).*epsilon_c))))),c_f.^( ...
-1).*(1+(-1).*beta).*delta.*(r+A.*theta.^(1+(-1).*B)+lambda).^(-1) ...
.*sigma,c_p.^(-1).*(1+(-1).*beta).*delta.*(r+A.*theta.^(1+(-1).*B) ...
+lambda).^(-1).*sigma];
return

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:

dF_depsilond

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

Dynamics

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-------------
BDG飽蠹閣 wechat
Enjoy it? Subscribe to my blog by scanning my public wechat account