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/
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
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.
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
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
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
@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
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.
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: