Domain error using nlsolve

I tried going through the Homotopy Continuation.jl website that introduced me to the basics of using the package. However the package seems to solve only polynomial systems. The non linear system of equations that my program uses has logarithmic functions too, which is why it is sensitive to the initial guess in the first place.
https://www.juliahomotopycontinuation.org/

Here’s the website I used for it

Thank you sir for your suggestion!

Could you please tell me how can Optim.jl be used to solve non-linear equation systems?
After going through the documentation of Optim.jl I find that it is used for finding the functional minimum and the corresponding points of minima
Thanks

Could you further elaborate on this sir? And could you please tell me more about pseudo arc continuation?

You will find many tutorials searching for something like “numerical homotopy continuation”. I find Appendix D of Boyd, J. P. (2001). Chebyshev and fourier spectral methods particularly nice.

Thanks for the reference sir! I shall check it out.

I think minimizing the squared residuals is similar to root finding. I am also exploring this package to find root:).

Ah! you mean setting the objective function to the sum of the squares of the resuidual functions?

Yes, but perhaps it is unnecessary to sum up.

If you need to find roots for a set of equations, then setting the objective functions to squared residuals of each equations may works.

1 Like

You can use PseudoArcLengthContinuation.jl as @Tamas_Papp advised to you, it works for any nonlinear function, and it scales well (> millions of unknowns).

2 Likes

Sure I will try it. Seems a bit complicated for me, but the reading the documentation should be fine. Could you cite some tutorials/examples? Will the Julia documentation be sufficient
Thank you so much for your reply

When I try installing PALC package as directed in the documentation: I receive the following error:

`(@v1.4) pkg>  add https://github.com/rveltz/PseudoArcLengthContinuation.jl
    Cloning git-repo `https://github.com/rveltz/PseudoArcLengthContinuation.jl`
   Updating git-repo `https://github.com/rveltz/PseudoArcLengthContinuation.jl`
  Resolving package versions...
ERROR: Unsatisfiable requirements detected for package RecipesBase [3cdcf5f2]:
 RecipesBase [3cdcf5f2] log:
 ├─possible versions are: [0.4.0, 0.5.0, 0.6.0, 0.7.0, 0.8.0] or uninstalled
 └─restricted to versions 1 by PseudoArcLengthContinuation [e5973c80] — no versions left
   └─PseudoArcLengthContinuation [e5973c80] log:
     ├─possible versions are: 0.1.0 or uninstalled
     └─PseudoArcLengthContinuation [e5973c80] is fixed to version 0.1.0`

I am running Julia version 1.4.

Interesting… I am following this announcement where it is asked to lower bound RecipesBase.

If you had RecipesBase installed, please update to 1.0

1 Like

Sorry for the late reply. I have updated RecipiesBase. And the installation was fine. Precompilation also smooth. Could I know the procedure using this package to find solutions to a system of non-linear equations?
Thank you so much

1 Like

So let’s say I have f,g,h as three equations whose solution I need to find. Then the objective function can be set to f^2+g^2+h^2 as the functional minimum is at f=g=h=0.

I tried writing out a simple code using the package. I used the newton solver as given in the first example.

using PseudoArcLengthContinuation
s_c = 5;#s/R_A2#reduced entropy
h_c = 0.75;#h/(R_A2*Θ_d)#reduced enthalpy
function f(x)
    F = similar(x)
    F[1] = h_c-((4+x[1])*x[2]+x[1])
    F[2] =  exp(x[1]-s_c) - (x[1]^(2*x[1]))*((1-x[1])^(1-x[1]))*((1/10^x[3])^(1+x[1]))/(x[2]^3)
    k = ((10^x[3]))*exp(-1/x[2])
    F[3] = x[1] - 0.5*(sqrt(k^2+4k)-k)
    return F
end 
x0 = [0.5, 0.05,6]
optnewton = NewtonPar(verbose=true)
out =@time newton(x->f(x),x0,optnewton)

I still get the same error

ERROR: LoadError: DomainError with -0.4561980584151881:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
Stacktrace:
 [1] throw_exp_domainerror(::Float64) at .\math.jl:37
 [2] ^ at .\math.jl:872 [inlined]
 [3] f(::Array{Float64,1}) at .\untitled-66927968220ce9667c9cc3493a5d6ab3:17
 [4] #19 at .\untitled-66927968220ce9667c9cc3493a5d6ab3:24 [inlined]
 [5] newton(::var"#19#20", ::PseudoArcLengthContinuation.var"#66#67"{var"#19#20"}, ::Array{Float64,1}, ::NewtonPar{Float64,DefaultLS,DefaultEig{typeof(real)}}; normN::typeof(norm), callback::Function, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\PseudoArcLengthContinuation\XY9lw\src\Newton.jl:103
 [6] newton(::Function, ::Function, ::Array{Float64,1}, ::NewtonPar{Float64,DefaultLS,DefaultEig{typeof(real)}}) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\PseudoArcLengthContinuation\XY9lw\src\Newton.jl:77
 [7] #newton#65 at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\PseudoArcLengthContinuation\XY9lw\src\Newton.jl:124 [inlined]
 [8] newton(::Function, ::Array{Float64,1}, ::NewtonPar{Float64,DefaultLS,DefaultEig{typeof(real)}}) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\PseudoArcLengthContinuation\XY9lw\src\Newton.jl:123
 [9] top-level scope at util.jl:175
in expression starting at untitled-66927968220ce9667c9cc3493a5d6ab3:24

I guess I need to go back and make my system more iterative-friendly!
Thanks for your help so far

It could be. However, if you have a parameter for which you know a solution, then you can use continuation to find the solutions for the parameters you are interested in

2 Likes

I do have a table of valid states that satisfy these equations. Let me try that too.

@PearsonCHEN I still end up with the same error using Optim.jl. I guess I need to modify my system of equations. Will close this for now :slight_smile:
Thanks for your help

My pleasure.

It could help if you slightly manipulate those equations. For example, you could divide the left-hand side by the right-hand side then minus 1. But as @Tamas_Papp said, you can not get rid of the problem of initial value easily when you use algorithms that find local optimum.

Good luck!

1 Like

Hello
I finally fixed this problem by running the function through a for loop, testing it at the points I was dealing with. Looks like modifying the initial guess based on the inputs suffices the problem. Thanks all of you to have shown me alternate options, which could come handy for my future endeavors :slight_smile:
Here’s the corrected function:

function mach(h,s,flag,a0,t0,r0)
    s_c = s/R_A2#reduced entropy
    h_c = h/(R_A2*Θ_d)#reduced enthalpy
    diff = h_c - (-2.0848e-18*s_c^9 + 1.8201e-15*s_c^8 + -6.6965e-13*s_c^7 +  1.3528e-10*s_c^6 +  -1.6406e-08*s_c^5 +  1.2298e-06*s_c^4 + -5.6915e-05*s_c^3 +  0.0016*s_c^2 + -0.0278*s_c + 0.4098)#dissociation line of 0.05 curve fitted
    #initilize values to avoid errors
    #consider a coarse lookup table for the initial guess
    if h_c<0.5
        a0 = 0.075;
        t0 =0.05;
        r0 =5;
    elseif h_c>1.125
        a0 = 0.8;
        t0 =0.1;
        r0 =6;
    elseif h_c >=0.5
        a0 =0.1;
        t0=0.07;
        r0 = 4;
    elseif h_c>=0.6 && s_c<7.35
        a0 =0.16;
        t0 =0.09;
        r0 =3;
    end
    #consider only dissociations above 5%
    if diff>=0
        function f!(F, x)
            F[1] =  h_c-((4+x[1])*x[2]+x[1])
            F[2] =  exp(x[1]-s_c) - (x[1]^(2*x[1]))*((1-x[1])^(1-x[1]))*((1/10^x[3])^(1+x[1]))/(x[2]^3)
            k = ((10^x[3]))*exp(-1/x[2])
            F[3] = x[1] - 0.5*(sqrt(k^2+4k)-k)
        end
        results = nlsolve(f!, [ a0,t0,r0])
        x = results.zero[1]#dissociation
        y = results.zero[2]#temperature ratio
        z = ρ_d/(10^results.zero[3])#density
        p = (1+x)*R_A2*z*y*Θ_d
        a = sqrt(R_A2*y*Θ_d*(x*(1-x^2)*(1+2*y) + (8+3*x-x^3)*y^2)/(x*(1-x) + 3*(2-x)*y^2))
        if flag==1
            return a
        else
            vec= [p,z,y*Θ_d,a,x]
        end
    else
        #println(" ")
        #println("Solver used: Perfect Gas Model")
        T = h/(4*R_A2)
        ρ = ρ_d*exp(-(s_c-3*log(T/Θ_d)))
        p = R_A2*T*ρ
        a = sqrt(4*R_A2*T/3)
        if flag==1
            return a
        else
            vec = [p,ρ,T,a,0]
        end
    end
end

So for those who might land up in this page cause they too faced a similar error:

  1. Try Optim.jl as suggested by @PearsonCHEN
  2. Try PALC as suggested by @rveltz and @Tamas_Papp
  3. If all else fails see where the system breaks down and if it rectifies with a simple change in initial guess.
    Thanks once again