Domain error using nlsolve

New to Julia
28

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