Solving a nonlinear system of equations in Julia


#1

I can write in function in Matlab in this way:

function res=resid(theta,alpha,beta);
RHS=[];
LHS=[];
RHS= theta-alpha;
LHS= theta*beta;
res = (LHS-RHS);

We set the parameters, call the function:

alpha=0.3;beta=0.95;
a01=[1.0;1.0];
th=fsolve(‘resid’,a01,[],alpha,beta)
This will return [6.0;6.0]. Does the option"[]" signal to fsolve that the input is a vector?

Anyway, how can I implement this in Julia using a NLsolve, Optim or JuMP? The original problem has more than 10 variables, so I would prefer a vector approach.

I can implement the function in Julia:

h! =function (theta)
RHS=[];
LHS=[];
RHS= theta-alpha;
LHS= theta*beta;
res= (LHS-RHS);
return res;
end
But simply using NLsolve:

a01 = [1.0;1.0];
res = nlsolve(h!,a01)
Returns:

MethodError: no method matching (::##17#18)(::Array{Float64,1}, ::Array{Float64,1})
Closest candidates are:
#17(::Any) at In[23]:3
If I alternatively use Optim, I get:

using Optim
optimize(h!, a01)
which returns:

MethodError: Cannot convert an object of type Array{Float64,1} to an object of type Float64
This may have arisen from a call to the constructor Float64(…),
since type constructors fall back to convert methods.


#2

A couple things to point out

  1. Semi colons are not necessary within functions
  2. In julia, when a function ends in an explamation point ! such as h!, that should signify that it mutates (changes) one of its inputs ( I say should since this is not enforced in the language but is a naming convention)
  3. When showing code on this forum, you should put your code in “backticks” (the key under tilda) for readability: this is enclosed in backticks

Now for the code, the error message you get from nlsolve has to do with point #2. NLsolve expects your function to take the residual vector and simply mutate it as opposed to returning a new vector. Here’s an example that works.

julia> using NLsolve

julia> function h!(theta, resvec, alpha, beta)
       resvec[1] = theta[1]*beta - (theta[1] - alpha)
       resvec[2] = theta[2]*beta - (theta[2] - alpha)
       end

julia> function f!(theta, resvec)
           h!(theta,resvec, 0.3, 0.95)
       end

julia> a01 = [1.0,1.0];

julia> nlsolve(f!,a01)
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.0,1.0]
 * Zero: [6.0,6.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 3
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 4
 * Jacobian Calls (df/dx): 4

I’m sure there are other ways to handle the alpha and beta parameters besides making another function.

Good reading materials:


https://docs.julialang.org/en/release-0.4/manual/noteworthy-differences/

ETA: And to clarify for your other example, Optim is used to optimize a function, this is different than solving it, the error you get is because you are returning a vector, while optimization is performed on a scalar (cost) quantity). Hope this helps!


#3

With Optim, from the error I see, I think you need to return just a number instead of an array. You didn’t specified exactly what you want to minimize, but I assume it’s some norm of resulting array. Try to return the norm https://docs.julialang.org/en/latest/stdlib/linalg/#Base.LinAlg.norm.


#4

You should cross-refer when double-posting: https://stackoverflow.com/questions/45545275/solving-a-nonlinear-system-of-equations-in-julia


#5

Hello!

Thanks everybody for helping!

Sorry for not cross-refering. Anyway, I posted this before getting to solution. Thanks ohsonice for your solution too.