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 
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
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:
- Try Optim.jl as suggested by @PearsonCHEN
- Try PALC as suggested by @rveltz and @Tamas_Papp
- If all else fails see where the system breaks down and if it rectifies with a simple change in initial guess.
Thanks once again