# Solve system of nonlinear equations

**URL:** <https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098>\
**Category:** New to Julia\
**Created:** [October 10, 2018, 3:34am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098 "2018-10-10T03:34:24Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![hirothreading](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hirothreading/32/5676_2.png) [@hirothreading](https://discourse.julialang.org/u/hirothreading)\
**Post date:** [October 10, 2018, 3:34am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/1 "2018-10-10T03:34:24Z")

</div>

H everyone

Apologies in advance for what may seem like a basic question, but I’m fairly new to julia. I’m looking to solve a system of nonlinear equations – currently I have 47 equations for 47 variables – and am looking for a solution more effective than MATLAB’s fsolve (unable to solve my problem) as I’m looking to find a steady state.

I think NLsolve would be most suitable for my needs, but I’m a bit unsure on some of the notation and syntax of julia.

Would the following format be appropriate for my problem:

> using NLsolve
> 
> function f!(F, x)  
> F[1] = …  
> F[2] = …  
> …  
> F[47] = …  
> end
> 
> nlsolve(f!, initial\_x, autodiff = :forward)

Where x is basically a vector of my endogenous variables, and initial\_x is a vector of starting guesses for the vector of endogenous variables.

Any help or advice would be much appreciated.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 10, 2018, 4:24am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/2 "2018-10-10T04:24:55Z")

</div>

Yes that looks right.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 10, 2018, 6:36am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/3 "2018-10-10T06:36:14Z")

</div>

> [@hirothreading](#):
>
> looking for a solution more effective than MATLAB’s fsolve (unable to solve my problem) as I’m looking to find a steady state

It is very unlikely that just switching to Julia will help you solve your problem, as the same general algorithms are used by pretty much all simple multivariate solvers. Usually the best approach is to look at the structure of your problem, and see if you can solve for blocks in variables (ie a block-triangular structure) independently of others, then break it up into simple ones. Good initial guesses also help a lot, which also need understanding of the problem.

Failing that, homotopy continuation methods may work, but ATM I am only aware of a Julia package for [polynomials](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl).

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [October 10, 2018, 7:53am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/4 "2018-10-10T07:53:23Z")

</div>

Yup, no general purpose continuation methods in Julia.

I think your (OP) best bet is to maybe show us the problem you’re trying to solve, because I think it’s very likely that NLsolve will fail it fsolve fails.

---

<div class="post-metadata">

**Author:** ![hirothreading](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hirothreading/32/5676_2.png) [@hirothreading](https://discourse.julialang.org/u/hirothreading)\
**Post date:** [October 10, 2018, 11:32am UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/5 "2018-10-10T11:32:50Z")

</div>

I have yet to code up the problem in julia, but I’ve just copied over the problem I’m trying to solve over from Atom. It’s basically a simple macroeconomic general equilibrium programming problem. Here’s the following system of equations:

```julia
            alpha = 1 - 0.623  
            n = 1.01^0.25 - 1  
            phiI = 132 
            kappa = 12 
            phiK = 1 
            delta = 0.025  
            omega = 1.023^(-0.25) 
            gamma = 1.1^(-0.25)  
            upsilon = 0.4  
            rho = -3  
            beta = .999 
            xi = 0.6  
            g = 0.17  
            b = 4  
            varrho = 0.4  
            phiPI = 2  
            phiR = 0.5  
            Rbar = (1/beta) - 1 
            phiNX = 1  
            ystar = 1 
            Rstar = (1/0.99 - 1) 
            phiS = 0.000742  
            fstarbar = -0.5  

function f!(F,x)
F[1]=x[2]^(alpha)*x[3]^(1-alpha) - x[1] 

F[2]=(x[5]/(1-alpha))^(1-alpha)*(x[6]/alpha)^alpha - x[4] 
     
F[3)=(1-alpha)*x[4]*x[1]/x[3] - x[5] 

F[4]=alpha*x[4]*x[1]/x[2] - x[6] 

F[5]=(1-x[4]-(phiI/2)*(x[7]-1)^2)*x[1] - x[10] 

F[6]=kappa*(x[4]-1)*x[1] + (1+n)*(x[7]/x[6])*phiI*(x[7]-1)*x[7]*x[1] - phiI*(x[7]-1)*x[7]*x[1]
     
F[7]=x[13]*x[12] - (1-phiK*((x[12]/x[12])-1)^2)*x[12] - x[11] 
     
F[8]=((1-delta)*x[2]+(1-(phiK*(x[12]/x[12])-1)^2)*x[12]) - x[2] 
     
F[9]=(x[7]/x[14])*(x[13]*(1-delta)+x[6]) - x[13] 
     
F[10]=(1-omega)+gamma*x[15] - (1+n)*x[15]
     
F[11]=(x[17]*x[18])^(-1/rho)*x[19]*(((1-upsilon)/upsilon)*(1/(xi*x[5])))^(1-upsilon) - x[16] 
     
F[12]=(x[17]*x[18])*(((x[14]*x[20])/x[7])+x[21]+x[22]-xi*x[25]) - x[19]
     
F[13]=x[23] + ((gamma*x[7])/x[14])*x[22] - x[22] 
     
F[14]=varrho*(x[5] - x[34]) - x[23] 
     
F[15]=x[15] - ((1-upsilon)/upsilon)*(1/(xi*x[5]))*x[19] - x[24] 
     
F[16]=xi*x[5]*x[24] + ((gamma*x[7])/x[14])*x[21] - x[21] 
     
F[17]=xi*x[25] + ((gamma*x[7])/x[14])*x[26] - x[26] 
     
F[18]=((x[14]*x[20])/x[7]) - x[19] + (1+omega)*(((x[14]*x[29])/x[7])+ x[5] - x[28]) - x[20] 
     
F[19]=1 - ((x[17]*x[18])/(x[17]*x[18]))*gamma*x[14]^(rho/(1-rho))*beta^(rho/(1-rho))*(1/x[7])^((upsilon*rho)/(1-rho))*(x[5]/x[5])^(((1-upsilon)*rho)/(1-rho)) - (x[17]*x[18])
     
F[20]=(x[18])^(-1/rho)*x[28]*(((1-upsilon)/upsilon)*(1/x[5]))^(1-upsilon) - x[27]
     
F[21]=x[18]*(((x[14]*x[29])/x[7])+ x[30] + x[31] - x[34]) - x[28]
     
F[22]=((x[7]*(x[32]-omega))/(x[32]*x[14]))*x[22] + ((x[7]*omega)/(x[32]*x[14]))*x[31] - x[31] 
     
F[23]=1 - ((1-upsilon)/upsilon)*(1/x[5])*x[28] - x[33]
     
F[24]=x[5]*x[33] + omega*((x[30]*x[7])/(x[32]*x[14])) + (1-omega)*(x[17])^((rho-1)/rho)*(1/xi)^(1-upsilon)*((x[21]*x[7])/(x[32]*x[14])) - x[30] 
     
F[25]=x[34] + omega*((x[35]*x[7])/(x[32]*x[14])) + (1-omega)*(x[17])^((rho-1)/rho)*(1/xi)^(1-upsilon)*((x[26]*x[7])/(x[32]*x[14])) - x[35] 
     
F[26]=omega*((x[14]*x[29])/x[7] + x[5]*x[33] - x[28]) - x[29] 
     
F[27]=1 - beta^(1/(1-rho))*(1/x[7])^((rho*upsilon)/(1-rho))*x[14]^(rho/(1-rho))*x[32]^(rho/(1-rho))*(x[5]/x[5])^(((upsilon-1)*rho)/(1-rho)) - x[18] 
     
F[28]=omega + (1-omega)*x[17]^(-(1-rho)/rho)*(1/xi)^(1-upsilon) - x[32] 
     
F[29]=x[19] + x[28] - x[36] 
     
F[30]=x[20]/x[37] - x[38] 
     
F[31]=omega*(1 - x[17]*x[18])*x[38]*((x[14]*A)/x[7]) - (1+n)*(x[38] - (1-omega))*x[37] 
     
F[32]=x[16] + x[27] - x[39] 
     
F[33]=(x[14]/x[7])*x[40] + x[41] + x[42] - x[43] - (1-n)*x[40] 
     
F[34]=x[34] + xi*x[25] - x[43] 
     
F[35]=g - x[41] 
     
F[36]=b - (1+n)*x[40] 
     
F[37]=varrho*PSI*((1-alpha)*x[4]*x[1] - x[43]) - x[42] 
     
F[38]=phiR*Rbar + (1-phiR)*x[14] + phiPI*x[7] - x[14] 

F[39]=(1+n)*(x[9]+x[11])/x[9] - x[14]/x[7] 
     
F[40]=(1+n)*(x[8]+x[10])/x[8] - x[14]/x[7] 
     
F[41]=x[33] + xi*x[24] - x[3] 
     
F[42]=x[36] + x[12] + x[41] + x[45] - x[1]*(1 - (phiI/2)*(x[7]-1)^2 - phiK*((x[12]/x[12])-1)^2) 
     
F[43]=x[44]^(phiNX)*ystar - x[45] 
     
F[44]=x[2] + x[40] + x[8] + x[44]*x[47] - x[37] 
     
F[45]=(1 + Rstar + x[46])*(x[44]/x[44]) - 1 - x[14] 
     
F[46]=(1 + Rstar + x[46])*x[44]*x[47] + x[45] - x[44]*x[47] 
     
F[47]=phiS*(exp(x[47]-fstarbar)-1) - x[46]
```

---

<div class="post-metadata">

**Author:** ![saschatimme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/saschatimme/32/10313_2.png) [@saschatimme](https://discourse.julialang.org/u/saschatimme)\
**Post date:** [October 10, 2018, 12:15pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/6 "2018-10-10T12:15:24Z")

</div>

I am a little bit confused, I see 23 things which look like parameters (`alpha`, `n`, `phiI`, …) and no variables (accesses to `x`). How do you get the 47 variables?

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [October 10, 2018, 12:54pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/7 "2018-10-10T12:54:31Z")

</div>

Variables are in all caps probably?

---

<div class="post-metadata">

**Author:** ![saschatimme](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/saschatimme/32/10313_2.png) [@saschatimme](https://discourse.julialang.org/u/saschatimme)\
**Post date:** [October 10, 2018, 1:02pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/8 "2018-10-10T13:02:04Z")

</div>

Oh that is it probably. So maybe one could formulate this as a polynomial system after some variable substitution and clearing of denominators.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 10, 2018, 1:03pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/9 "2018-10-10T13:03:19Z")

</div>

> [@hirothreading](#):
>
> I have yet to code up the problem in julia

Usually, thinking about the mathematical structure of the problem first (before programming) has large and immediate payoffs. This looks like Euler and technology equation residuals from a DSGE model, so maybe you can derive the steady state capital, then consumption and labor, etc.

---

<div class="post-metadata">

**Author:** ![hirothreading](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hirothreading/32/5676_2.png) [@hirothreading](https://discourse.julialang.org/u/hirothreading)\
**Post date:** [October 10, 2018, 1:10pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/10 "2018-10-10T13:10:02Z")

</div>

Apologies, the code was essentially pulled out of a MATLAB script. I’ll edit my post now.

---

<div class="post-metadata">

**Author:** ![hirothreading](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/hirothreading/32/5676_2.png) [@hirothreading](https://discourse.julialang.org/u/hirothreading)\
**Post date:** [October 10, 2018, 1:13pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/11 "2018-10-10T13:13:12Z")

</div>

Correct, it’s a DSGE model.

Worst case scenario, I may have to manually log linearise the model and put into Dynare, setting initial values equal to zero.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [October 10, 2018, 1:38pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/12 "2018-10-10T13:38:51Z")

</div>

As I said, it should be possible to simplify a lot by understanding the model. Brute-forcing your way to a steady state will be costly (it is actually the most costly part of solving large-ish DSGE models by perturbation methods; second is the linear part, the rest is pretty mechanical).

You can find examples on how to do this in most texts, eg the handbook chapter of Fernández-Villaverde et al. (Sorry if you already know this).

---

<div class="post-metadata">

**Author:** ![aaowens](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aaowens/32/12101_2.png) [@aaowens](https://discourse.julialang.org/u/aaowens)\
**Post date:** [October 10, 2018, 4:13pm UTC](https://discourse.julialang.org/t/solve-system-of-nonlinear-equations/16098/13 "2018-10-10T16:13:46Z")

</div>

One thing to consider is the range of the valid parameter space. If you have some variables, like variances, which have to be positive, solve for the log of that variable and exp it in the objective.

Related, if there are outputs like consumption and you have log utility, then the marginal utility 1/c only makes sense for positive consumption. You might try truncating consumption to 0 so you get an Inf (and fail on that particular guess), instead of crossing the discontinuity at 0, which might make the optimizer diverge.