Hi, thanks for writing and publishing the package!
In evaluating EquationsSolver.jl
package, I tried to solve a problem arising in balancing a chemical reaction with chemical equilibrium equations as following (a Textbook gives the solution as (x, y, z, w) = (0.271, 0.213, 1.849, 1.032)):
julia> @variables x y z w;
julia> eqs = [
2 ~ 2x + 2y + w,
5 ~ x + 2z + w,
0.5816 ~ (y * z^(1//2) / x) * inv(x + y + w + z)^(1//2),
0.9570 ~ (w * y^(1//2) / x) * inv(x + y + w + z)^(1//2)
];
julia> vars = Dict(x=>1, y=>1, z=>1, w=>1);
julia> pro = NLProblem(eqs,vars);
julia> res = solve(pro)
ERROR: InexactError: Int64(0.0013379310344827644)
...
I then figured I should have declared vars
with Float64
's and not Int
's, and ran:
julia> vars = Dict(x=>1.0, y=>1.0, z=>1.0, w=>1.0);
julia> pro = NLProblem(eqs,vars);
julia> res = solve(pro)
ERROR: DomainError with -0.18794749606644048:
Exponentiation yielding a complex result requires a complex argument.
...
I suspect the solving algorithm wanders in Complex
terrain before converging to the real solution (given in the Textbook), with very low imaginary components. I then declared vars
as Complex
, by running:
julia> vars = Dict(x=>1.0+0.0im, y=>1.0+0.0im, z=>1.0+0.0im, w=>1.0+0.0im);
julia> pro = NLProblem(eqs,vars);
julia> res = solve(pro)
ERROR: InexactError: Float64(0.0 + 0.08218679783487991im)
Stacktrace:
[1] Real
@ ./complex.jl:44 [inlined]
[2] convert
@ ./number.jl:7 [inlined]
[3] setindex!(A::Matrix{Float64}, x::ComplexF64, i1::Int64)
@ Base ./array.jl:903
[4] _unsafe_copyto!(dest::Matrix{Float64}, doffs::Int64, src::Matrix{Number}, soffs::Int64, n::Int64)
@ Base ./array.jl:253
[5] unsafe_copyto!
@ ./array.jl:307 [inlined]
[6] _copyto_impl!
@ ./array.jl:331 [inlined]
[7] copyto!
@ ./array.jl:317 [inlined]
[8] copyto!
@ ./array.jl:343 [inlined]
[9] copyto_axcheck!(dest::Matrix{Float64}, src::Matrix{Number})
@ Base ./abstractarray.jl:1104
[10] AbstractArray
@ ./array.jl:564 [inlined]
[11] AbstractArray
@ ./boot.jl:484 [inlined]
[12] convert
@ ./abstractarray.jl:17 [inlined]
[13] inv(A::Matrix{Number})
@ LinearAlgebra /usr/share/julia/stdlib/v1.7/LinearAlgebra/src/dense.jl:870
[14] one_cal(array::Vector{ComplexF64}, dict::Dict{Symbolics.Num, ComplexF64}, f::Vector{Symbolics.Num}, ja::Matrix{Symbolics.Num})
@ EquationsSolver ~/.julia/packages/EquationsSolver/sjbmW/src/solver.jl:13
[15] _solve(dict::Dict{Symbolics.Num, ComplexF64}, f::Vector{Symbolics.Num}, ja::Matrix{Symbolics.Num}, maxiter::Int64, abstol::Float64, vars::Vector{Symbolics.Num})
@ EquationsSolver ~/.julia/packages/EquationsSolver/sjbmW/src/solver.jl:33
[16] solve(problem::EquationsSolver.NonlinearProblem)
@ EquationsSolver ~/.julia/packages/EquationsSolver/sjbmW/src/solver.jl:54
[17] top-level scope
@ REPL[35]:1
Was the package written with Complex
variables in mind?
FINALLY then, I tried much closer initial values (that I wouldn't be able to guess if I didn't have the solution from the Textbook):
julia> vars = Dict(x=>0.3, y=>0.2, z=>1.8, w=>1.0);
julia> pro = NLProblem(eqs,vars);
julia> res = solve(pro)
Dict{Symbolics.Num, Float64} with 4 entries:
y => 0.212802
x => 0.271204
w => 1.03199
z => 1.8484 # The Textbook values at last!
Are there ways to control the solver as to avoid what it seems to be a diverging pattern of the cases above? (Such as relaxation factors)?