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.
You can use PseudoArcLengthContinuation.jl as @Tamas_Papp advised to you, it works for any nonlinear function, and it scales well (> millions of unknowns).
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
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
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!
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: