Domain error using nlsolve

I am trying to solve non linear equation systems in gas dynamics using nlsolve, but the solver seems to breakdown when I call it at h = 0.1R_A2Theta_D (see code below)

using NLsolve
using DifferentialEquations
## About this code module
# Julia Program to solve Prandtl-Meyer flow considering high temperature effects
## GAS PROPERTIES
#comment out the appropriate section/lines to select
##Oxygen
R_A2 = 8.314*1000/32; #the molar gas constant for Oxygen
Θ_d = 59500; #dissociation temperature in Kelvin
ρ_d = 150000; #dissociation density
#flag = 1; #to display on the plot
##Nitrogen
#R_A2 = 8.314*1000/28; #the molar gas constant for Oxygen
#Θ_d = 113000; #dissociation temperature in Kelvin
#ρ_d = 130000; #dissociation density
#flag =2; #to display on the plot
##SOLUTION
## STATE EQUATION for high temperature gas dynamics in-terms of pressure and temperature input
#Intended to calculate inlet properties
function state(P,T)
#x[1]-->α_*, x[2]-->log10(ρ_d/ρ)
    function g!(G,x)
        G[1] = P - (1+x[1])*R_A2*T*ρ_d/10^x[2]
        k = (10^(x[2]))*exp(-Θ_d/T)
        G[2] = x[1] - 0.5*(sqrt(k^2 + 4k) -k)
    end
    res = nlsolve(g!,[0.5,6])
    x = res.zero[1];#dissociation
    y = res.zero[2];#density ratio
    s = R_A2*((3*log(T/Θ_d) + (x*(1-2*log(x)))-((1-x)*log(1-x)) - (1+x)*log(1/10^y)))
    a = sqrt(R_A2*(T/Θ_d)*Θ_d*(x*(1-x^2)*(1+2*(T/Θ_d)) + (8+3*x-x^2)*(T/Θ_d)^2)/(x*(1-x) + 3*(2-x)*(T/Θ_d)^2))
    h = R_A2*((4+x)*T + x*Θ_d)
    #comment out the following two lines to get absolute enthalpies
    #h = h/(R_A2*Θ_d)
    #s = s/R_A2
    vec = [s,a,h,T/Θ_d,log10(R_A2*ρ_d*Θ_d/P),y,x]
    return vec
end

##Mach function calculator for High-Temperature Gas Dynamics
#NON-LINEAR EQUATION SYSTEM FOR STATE
#x[1]-->α_*, x[2]-->T/Θ_d, x[3]-->log10(ρ_d/ρ)
function mach(h ,s)
    s_c = s/R_A2#reduced entropy
    h_c = h/(R_A2*Θ_d)#reduced enthalpy
    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!, [ 0.5, 0.01,4.5])
    x = results.zero[1]#dissociation
    y = results.zero[2]#temperature ratio
    z = ρ_d/(10^results.zero[3])#density
    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))
    vec = [x,y,z,a]
    return vec
end
## Main program for Prandtl-Meyer Expansion
##Upstream properties
p_a = 121590;#upstream pressure in pascal
T_a = 6140;#upstream temperature in kelvin
M_a = 1;#upstream Mach number
state_info = state(p_a,T_a);
s_a = state_info[1];#inlet entropy
u_a = M_a*state_info[2];#inlet velocity
h_a = state_info[3];#inlet enthalpy
h_t_a = h_a + 0.5*u_a^2;#inlet stagnaion enthalpy
##Prandtl-Meyer Solver
res = mach(0.1*R_A2*Θ_d,s_a)
println("  ");
println(res);

I get the following error:

LoadError: DomainError with -6.055454134369954e-6:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
in expression starting at C:\Users\SURYA\Desktop\IIST Internship\Programs\Julia programs\MachNumber.jl:99
throw_exp_domainerror(::Float64) at math.jl:37
^ at math.jl:872 [inlined]
(::var"#f!#66"{Float64,Float64})(::Array{Float64,1}, ::Array{Float64,1}) at MachNumber.jl:62
finite_difference_jacobian!(::Array{Float64,2}, ::var"#f!#66"{Float64,Float64}, ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}, ::Nothing; relstep::Float64, absstep::Float64, colorvec::UnitRange{Int64}, sparsity::Nothing, dir::Bool) at jacobians.jl:370
finite_difference_jacobian!(::Array{Float64,2}, ::Function, ::Array{Float64,1}, ::FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}, ::Nothing) at jacobians.jl:295
finite_difference_jacobian! at jacobians.jl:295 [inlined]
fj_finitediff! at oncedifferentiable.jl:139 [inlined]
(::NLSolversBase.var"#j_finitediff!#22"{Array{Float64,1},NLSolversBase.var"#fj_finitediff!#21"{var"#f!#66"{Float64,Float64},FiniteDiff.JacobianCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}}})(::Array{Float64,2}, ::Array{Float64,1}) at oncedifferentiable.jl:144
jacobian!!(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,2}, ::Array{Float64,1}) at interface.jl:142
jacobian!! at interface.jl:140 [inlined]
jacobian! at interface.jl:135 [inlined]
trust_region_(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Array{Float64,1}}) at trust_region.jl:181
trust_region at trust_region.jl:229 [inlined]
trust_region at trust_region.jl:229 [inlined]
nlsolve(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64) at nlsolve.jl:26
nlsolve at nlsolve.jl:18 [inlined]
nlsolve(::Function, ::Array{Float64,1}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at nlsolve.jl:52
nlsolve(::Function, ::Array{Float64,1}) at nlsolve.jl:46
mach(::Float64, ::Float64) at MachNumber.jl:66
top-level scope at MachNumber.jl:99

It looks like a problem of initial guess.

Can you try this?

    results = nlsolve(f!, [0.5, 0.01,10], show_trace=true, ftol = 1e-6)

This guess gives the following result:

Iter     f(x) inf-norm    Step 2-norm 
------   --------------   --------------
     0     4.999986e-01              NaN
     1     5.561852e-03     3.302127e+00
     2     2.065245e-05     3.597072e-01
     3     1.252077e-04     0.000000e+00
     4     1.252077e-04     0.000000e+00
     5     1.252077e-04     0.000000e+00
     6     4.893049e-05     0.000000e+00
     7     1.274549e-05     1.329266e-01
     8     4.168710e-05     0.000000e+00
     9     9.916628e-06     1.281029e-01
    10     8.574259e-06     1.211173e-01
    11     7.406114e-06     1.149321e-01
    12     6.432701e-06     1.094965e-01
    13     5.611286e-06     1.046729e-01
    14     4.912477e-06     1.003605e-01
    15     4.313567e-06     9.647944e-02
    16     3.796859e-06     9.296623e-02
    17     3.348377e-06     8.976929e-02
    18     2.956960e-06     8.684643e-02
    19     2.613608e-06     8.416273e-02
    20     2.311000e-06     8.168906e-02
    21     2.043126e-06     7.940086e-02
    22     1.805006e-06     7.727734e-02
    23     1.592453e-06     7.530072e-02
    24     1.401843e-06     7.345571e-02
    25     1.229684e-06     7.172908e-02
    26     4.363853e-07     4.480841e-02
  
[0.010112658920948703, 0.022415164678805527, 6.121796865163041e-11, 623.5857844645709]

Hello
Thanks for your reply

If I change the initial guess, it works fine for this case, but throws am error for other pairs of h and s. After experimenting with a lot of pairs of h and s, one can conclude the system is sensitive to the initial guess. I am looking for some globally convergent methods like Homotopy Continuation but in Julia, it works only for polynomial systems. Would appreciate if you have any other way around it.

Thanks once again

with this initial guess, the system blows up at
some other (h,s). Look at the following case for example:

using NLsolve
using Plots
#Oxygen
R_A2 = 8.314*1000/32; #the molar gas constant for Oxygen
Θ_d = 59500; #dissociation temperature in Kelvin
ρ_d = 150000; #dissociation density
#Nitrogen
#R_A2 = 8.314*1000/28; #the molar gas constant for Nitrogen
#Θ_d = 113000; #dissociation temperature in Kelvin
#ρ_d = 130000; #dissociation density
#γ = 1.33; #  for calculation of initial guess
#c_p = 3*R_A2;; #  for calculation og initial guess
#state function
function mach(h,s,flag,a0,t0,r0)
    s_c = s/R_A2#reduced entropy
    h_c = h/(R_A2*Θ_d)#reduced enthalpy
    #initialize guess based on location in chart
    
    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
    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], show_trace=true, ftol = 1e-6)
        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= [log10((R_A2*Θ_d*ρ_d)/p),log10(ρ_d/z),y,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
h = 1.35*R_A2*Θ_d;
s = 5*R_A2
println(" h = ", h/(R_A2*Θ_d), " s = ", 5)
res = mach(h,s,2,0.5, 0.01,10);  
println(res)

This is a revised function of the same, where I have restricted the use of the nonlinear solver for values (h,s) where I am certain it will diverge using a simple curvefit. Despite, other points seem to diverge because of the initial guess as you have pointed out rightly. When I run the code above I get the following output:

h = 1.351.3500000000000003 s = 5
Iter     f(x) inf-norm    Step 2-norm 
------   --------------   --------------
     0     8.050000e-01              NaN
     1              Inf     8.673617e-18
     2              Inf     5.551115e-17
     3              Inf     5.551115e-17
     4              Inf     0.000000e+00
     5              Inf     0.000000e+00
     6              Inf     0.000000e+00
     7              Inf     0.000000e+00
     8              Inf     0.000000e+00
     9              Inf     0.000000e+00
    10              Inf     0.000000e+00
    11              Inf     5.551115e-17
    12              Inf     0.000000e+00
    13              Inf     5.551115e-17
    14              Inf     5.551115e-17
    15              Inf     5.551115e-17
    16              Inf     0.000000e+00
    17    1.584226e+168     5.551115e-17
    18     2.366336e+79     0.000000e+00
    19     9.146412e+34     0.000000e+00
    20     5.686554e+12     5.551115e-17
    21     4.482740e+01     5.551115e-17
    22     8.049996e-01     3.701072e+00
    23     8.012912e-01     9.742420e-02
    24     7.938775e-01     1.721280e-01
    25     7.790625e-01     2.801262e-01
    26     7.494888e-01     4.136966e-01
    27     6.906168e-01     5.587927e-01
    28     5.744490e-01     1.543148e-01
    29     4.568162e-01     6.400819e+00
    30     2.610001e-01     2.041366e-01
    31     1.831104e-02     2.692114e-01
    32              Inf     1.563194e-13
    33              Inf     1.387779e-17
    34              Inf     0.000000e+00
    35              Inf     0.000000e+00
    36              Inf     0.000000e+00
    37              Inf     0.000000e+00
ERROR: LoadError: DomainError with -0.00028739291806911815:
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] (::var"#f!#323"{Float64,Float64})(::Array{Float64,1}, ::Array{Float64,1}) at .\untitled-594fe38f86d51795b3993c1e9cdc3dc3:23
 [4] value!!(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLSolversBase\mGaJg\src\interface.jl:166
 [5] value!! at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLSolversBase\mGaJg\src\interface.jl:164 [inlined]
 [6] value! at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLSolversBase\mGaJg\src\interface.jl:28 [inlined]
 [7] trust_region_(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Array{Float64,1}}) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLsolve\ZBTu4\src\solvers\trust_region.jl:171
 [8] trust_region at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLsolve\ZBTu4\src\solvers\trust_region.jl:229 [inlined] (repeats 2 times)
 [9] nlsolve(::OnceDifferentiable{Array{Float64,1},Array{Float64,2},Array{Float64,1}}, ::Array{Float64,1}; method::Symbol, xtol::Float64, ftol::Float64, iterations::Int64, store_trace::Bool, show_trace::Bool, extended_trace::Bool, linesearch::Static, linsolve::NLsolve.var"#27#29", factor::Float64, autoscale::Bool, m::Int64, beta::Int64, aa_start::Int64, droptol::Float64) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLsolve\ZBTu4\src\nlsolve\nlsolve.jl:26
 [10] nlsolve(::Function, ::Array{Float64,1}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol},NamedTuple{(:show_trace, :ftol),Tuple{Bool,Float64}}}) at C:\Users\SURYA\.juliapro\JuliaPro_v1.4.0-1\packages\NLsolve\ZBTu4\src\nlsolve\nlsolve.jl:52
 [11] mach(::Float64, ::Float64, ::Int64, ::Float64, ::Float64, ::Int64) at .\untitled-594fe38f86d51795b3993c1e9cdc3dc3:27
 [12] top-level scope at untitled-594fe38f86d51795b3993c1e9cdc3dc3:55
in expression starting at untitled-594fe38f86d51795b3993c1e9cdc3dc3:55

There are no general black-box methods for systems of nonlinear equations. You have to invest a bit into understanding your problem. If you can find a well-behaved special case,

can be used for a homotopy continuation.

Thanks for the suggestion.
Yes you are right, I have just thought about putting some if conditions on the input and modify the initial guess values. Using homotopy continuation seems far fetched as of now!

Besides, you can try different algorithms in NLsolve.jl, different packages like Optim.jl, with auto differentiation and line search.

There is the “naive homotopy” approach of just using solutions from a succession of parametrized problems as the initial guess of the next one (ie w/o derivatives). They are easy to program but do not scale well. Nevertheless, they can be better than nothing.

Homotopy continuation is not that scary and worth investing in if the problem is differentiable — Julia has very accessible AD solutions so it is not hard to implement.

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