# Calling Mathematica into Julia to symbolically solve a system of non-linear equations

I’m currently using Mathematica to symbolically solve a system of non-linear ODEs (10 equations with 8 variables) as I believe this is not yet possible using Symbolics.jl. Is anyone able to comment on whether it would be possible to call Mathematica into Julia using MathLink.jl to solve the below system. The documentation provided a simple example that I wasn’t able to extrapolate to my Mathematica code (below) successfully. Thanks very much.

The system is governed by 14 parameters:

``````kb1 = 0.01;
K1  = 10.0;
d1  = 1.0;
l11 = 1.0;
m11 = 0.01;
l12 = 1.0;
m12 = 5.0;
kb2 = 0.01;
K2  = 30.0;
d2  = 1.0;
l22 = 1.0;
m22 = 0.01;
l21 = 1.0;
m21 = 5.0;
``````

and is solved in the following way in Mathematica:

``````Solve[{m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0,
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0,
l11*M1*X10 - m11*X11 == 0,
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0,
l12*M2*X10 - m12*X12 == 0,
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0,
l22*M2*X20 - m22*X22 == 0, l21*M1*X20 - m21*X21 == 0,
X10 + X12 + X11 == 1,
X20 + X21 + X22 == 1
&& M1 >= 0
&& M2 >= 0} ,
{X10, M1, X11, M2, X12, X20, X22, X21}, Reals]
``````

generating the below fixed points for variables M1 and M2:

``````{{X10 -> 0.000999999, M1 -> 9.93007, X11 -> 0.993006, M2 -> 29.9701,
X12 -> 0.00599402, X20 -> 0.000333333, X22 -> 0.999005,
X21 -> 0.000662005}}
``````

How could I obtain this result in Julia?

1 Like

Seems to me you are simply solving a system of linear equations here. You can use `JuMP` for that.

It’s nonlinear, because there are terms like `M1*X10` that are products of unknowns. However, I don’t see how it is an “ODE” as claimed by @MAC … I guess it comes from the fixed points of an ODE? — it looks like a system of polynomial equations, which can be solved very efficiently by the HomotopyContinuation.jl package if you want a numerical answer (which is what you are getting from Mathematica anyway, it seems … in what way are you solving it “symbolically”?).

4 Likes

This is over-determined, so may not have a solution, but clearly it does since you provide one.

Here is how you might do it with JuMP and the Gurobi solver:

JuMP Example
``````using JuMP
using Gurobi

## Parameters:
kb1 = 0.01;
K1  = 10.0;
d1  = 1.0;
l11 = 1.0;
m11 = 0.01;
l12 = 1.0;
m12 = 5.0;
kb2 = 0.01;
K2  = 30.0;
d2  = 1.0;
l22 = 1.0;
m22 = 0.01;
l21 = 1.0;
m21 = 5.0;

## Model def.:
model = JuMP.Model(Gurobi.Optimizer)
JuMP.set_optimizer_attribute(model, "NonConvex", 2)

## Variables:
V = @variables(model, begin
X10
M1
X11
M2
X12
X20
X22
X21
end)

C = @constraints(model, begin
m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0
l11*M1*X10 - m11*X11 == 0
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0
l12*M2*X10 - m12*X12 == 0
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0
l22*M2*X20 - m22*X22 == 0
l21*M1*X20 - m21*X21 == 0
X10 + X12 + X11 == 1
X20 + X21 + X22 == 1
M1 >= 0
M2 >= 0
end)

optimize!(model)

result_dict = Dict(V .=> value.(V))
primal_feasibility_report(model, result_dict)

# Provided solution:
feas_pt = Dict(X10 => 0.000999999, M1 => 9.93007, X11 => 0.993006, M2 => 29.9701,
X12 => 0.00599402, X20 => 0.000333333, X22 => 0.999005,
X21 => 0.000662005)

primal_feasibility_report(model, feas_pt)
``````
4 Likes

There are 8 equality constraints in 8 variables. Whoops, I miscounted — 10 equality constraints in 8 variables. There are two inequalities which filter out some (in this case all but one) of the roots.

With HomotopyContinuation.jl I can reproduce the Mathematica solution:

``````julia> using HomotopyContinuation

julia> kb1 = 0.01;
julia> K1  = 10.0;
julia> d1  = 1.0;
julia> l11 = 1.0;
julia> m11 = 0.01;
julia> l12 = 1.0;
julia> m12 = 5.0;
julia> kb2 = 0.01;
julia> K2  = 30.0;
julia> d2  = 1.0;
julia> l22 = 1.0;
julia> m22 = 0.01;
julia> l21 = 1.0;
julia> m21 = 5.0;

julia> @var X10 M1 X11 M2 X12 X20 X22 X21
(X10, M1, X11, M2, X12, X20, X22, X21)

julia> F = System([m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10,
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20,
l11*M1*X10 - m11*X11,
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20,
l12*M2*X10 - m12*X12,
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20,
l22*M2*X20 - m22*X22 == 0, l21*M1*X20 - m21*X21,
X10 + X12 + X11 - 1,
X20 + X21 + X22 - 1])
System of length 10
8 variables: M1, M2, X10, X11, X12, X20, X21, X22

0.01*X11 + 5.0*X12 - 1.0*M1*X10 - 1.0*M2*X10
-1.0*M1 + 0.01*X10 + 10.01*X11 + 5.0*X21 - 1.0*M1*X10 - 1.0*M1*X20
-0.01*X11 + 1.0*M1*X10
-1.0*M2 + 5.0*X12 + 0.01*X20 + 30.01*X22 - 1.0*M2*X10 - 1.0*M2*X20
-5.0*X12 + 1.0*M2*X10
5.0*X21 + 0.01*X22 - 1.0*M1*X20 - 1.0*M2*X20
0
-5.0*X21 + 1.0*M1*X20
-1 + X10 + X11 + X12
-1 + X20 + X21 + X22

julia> result = solve(F)
Tracking 6 paths... 100%|███████████████████████████████| Time: 0:00:16
# paths tracked:                  6
# non-singular solutions (real):  4 (4)
# singular endpoints (real):      0 (0)
# total solutions (real):         4 (4)
Result with 4 solutions
=======================
• 6 paths tracked
• 4 non-singular solutions (4 real)
• random_seed: 0x86ff6d6b
• start_system: :polyhedral

julia> real_solutions(result) # extract the real roots
4-element Vector{Vector{Float64}}:
[-1.0070462962057125e-5, 29.990003354585465, 0.142918523818753, -0.00014392557011470473, 0.8572254017513615, 0.00033333329628380104, -6.713646621069738e-10, 0.9996666673750833]
[9.930069784028507, 29.97014319708603, 0.0009999989929586964, 0.9930059784038578, 0.005994022603182795, 0.0003333332962592715, 0.0006620045786384862, 0.9990046621251022]
[9.9900100166733, -3.3366674193768253e-6, 0.0009999989990009555, 0.9990000016683313, -6.67332813376323e-10, 0.33359260855930595, 0.6665187001991278, -0.00011130875843388742]
[-1.0009999973788885e-5, -3.3344444459124908e-6, 1.001002671228785, -0.0010020036712260715, -6.675575594819791e-7, 1.0003355590062681, -2.002671784362474e-6, -0.00033355633448375795]

julia>  filter(x -> x[1] >= 0 && x[2] >= 0, real_solutions(result)) # extract the M1 >= 0, M2 >= solutions
1-element Vector{Vector{Float64}}:
[9.930069784028507, 29.97014319708603, 0.0009999989929586964, 0.9930059784038578, 0.005994022603182795, 0.0003333332962592715, 0.0006620045786384862, 0.9990046621251022]
``````

Note that the final line is the solution for `M1, M2, X10, X11, X12, X20, X21, X22` (`== F.variables`), and seems to match what you are reporting from Mathematica.

8 Likes

?

HomotopyContinuation.jl is great too.

Whoops, right, I miscounted commas when scanning the equations by eye. Luckily, Julia’s parser is better at counting than me.

1 Like

I’m curious how long Gurobi took to solve this? HomotopyContinuation took 16 seconds (though it was finding all the roots, and I only filtered them a posteriori). Looks like this is mostly compilation time; solving a second time is about 0.02 seconds.

Initial:
`7.042219 seconds (19.25 M allocations: 1.061 GiB, 5.44% gc time, 99.02% compilation time)`
After compiliation:
`0.058017 seconds (1.80 k allocations: 96.125 KiB)`

This is finding a (slightly better polished) point that matches the given one.

Oh, so it’s not really comparable — it’s taking an approximate solution to start with, whereas HomotopyContinuation is finding the solution from scratch?

(@MAC, I wonder how HomotopyContinuation’s performance compares to Mathematica here?)

Sorry if that wasn’t clear, it is starting from scratch, just the result satisfies the constraints to a smaller tolerance than the provided one.

2 Likes

Since this is nonconvex, is Gurobi guaranteed to find the solutions if they exist?

Yes, this setting in the newer versions of Gurobi guarantees global lower bounds and optimality for bilinear problems. That said, this is a feasbility problem. I also ran it asking it to find all available solutions, and it returns this point only. It becomes much harder if we drop the lower bounds `M1 >= 0, M2 >= 0` and filter a posteriori. I think HomotopyContinuation is more impressive here in finding multiple solutions.

3 Likes

@jd-foster and @stevengj: Thanks so much for all these great responses - it’s my first time posting here and it has been very helpful.

In the Mathematica code the system is solved exactly and the results are subsequently numericized - sorry for the confusion. A numerical solution will still work for me so these examples are excellent. The benchmarking is helpful too as I will run the system over a very large parameter set.

@stevengj in Mathematica this example takes on average 0.0097s to solve.

``````RepeatedTiming[
Solve[{m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0,
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0,
l11*M1*X10 - m11*X11 == 0,
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0,
l12*M2*X10 - m12*X12 == 0,
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0,
l22*M2*X20 - m22*X22 == 0,
l21*M1*X20 - m21*X21 == 0,
X10 + X12 + X11 == 1,
X20 + X21 + X22 == 1
&& M1 >= 0 && M2 >= 0} ,
{X10, M1, X11, M2, X12, X20, X22, X21}, Reals], 30]
``````
``````{0.0097, {{X10 -> 0.000999999, M1 -> 9.93007, X11 -> 0.993006, M2 -> 29.9701, X12 -> 0.00599402, X20 -> 0.000333333, X22 -> 0.999005, X21 -> 0.000662005}}}
``````

If you are varying parameters and computing roots over and over, the problem becomes much easier, because then you just need to track how solutions evolve. (This is called a “homotopy” or “numerical continuation”. HomotopyContinuation has built-in methods for this.)

Even without taking advantage of this, in HomotopyContinuations it looks like much of the 16 seconds I quote above was compilation time (or some other initialization? I don’t know how much it is caching things internally). Solving again is about 0.02 seconds.

3 Likes

As an aside: I ran a different parameter set that should produce 3 fixed points (2 stable and 1 unstable) as determined by Mathematica. Using JuMP with Gurobi it was only able to identify 1 fixed point while HomotopyContinuation was able to find all 3 (unless I have implemented something incorrectly).

``````#Parameters
kb1 =  0.01;
K1  =  50.0;
d1  =  1.0;
l11 =  1.0;
m11 =  10.0;
l12 =  1.0;
m12 =  1.0;
kb2 =  0.01;
K2  =  20.0;
d2  =  1.0;
l22 =  5.0;
m22 =  1.0;
l21 =  10.0;
m21 =  1.0;
``````
``````using HomotopyContinuation

@var X10 M1 X11 M2 X12 X20 X22 X21

F = System([m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10),
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20,
l11*M1*X10 - m11*X11,
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20,
l12*M2*X10 - m12*X12,
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20,
l22*M2*X20 - m22*X22,
l21*M1*X20 - m21*X21,
X10 + X12 + X11 - 1,
X20 + X21 + X22 - 1])

result = HomotopyContinuation.solve(F)

#Finds all 3 fixed points (M1=>40, M2=>0); (M1=>8.32, M2=>3.17); (M1=>0, M2=>19.8)
all_fps = filter(x -> x[1] >= 0 && x[2] >= 0, real_solutions(result))
``````
``````using JuMP
using Gurobi

model = JuMP.Model(Gurobi.Optimizer)
JuMP.set_optimizer_attribute(model, "NonConvex", 2)

V = @variables(model, begin
X10
M1 >= 0
X11
M2 >= 0
X12
X20
X22
X21
end)

C = @constraints(model, begin
m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0
l11*M1*X10 - m11*X11 == 0
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0
l12*M2*X10 - m12*X12 == 0
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0
l22*M2*X20 - m22*X22 == 0
l21*M1*X20 - m21*X21 == 0
X10 + X12 + X11 == 1
X20 + X21 + X22 == 1
end)

optimize!(model)

#Seems to find only 1 fixed point M1=>40, M2=>0
result_dict = Dict(V .=> value.(V))
``````

Hi Megan,

``````JuMP.set_optimizer_attribute(model, "PoolSearchMode", 2)
JuMP.set_optimizer_attribute(model, "PoolSolutions", 1000)
``````

where the 1000 in `"PoolSolutions"` is some arbitrary but large enough number for the number of solutions you expect to capture.
The default in Gurobi is to return just one optimal solution, and since the objective is just 0 here, it returns one feasible point.
With the above, I get 14 solutions, but these are not all different up to small tolerances.

The 'unique' solutions are:
``````  X10 => 0.19999
X22 => 4.14194e-7
X11 => 0.800003
X20 => 0.00249363
M2  => 3.32202e-5
X21 => 0.997506
X12 => 6.6437e-6
M1  => 40.0022
``````
`````` X10 => 0.0480795
X22 => 0.989937
X11 => 3.04079e-6
X20 => 0.00999995
M2  => 19.7988
X21 => 6.32807e-5
X12 => 0.951917
M1  => 0.00063283
``````
``````  X10 => 0.199952
X22 => 0.158478
X11 => 0.16627
X20 => 0.00999969
M2  => 3.16965
X21 => 0.831523
X12 => 0.633778
M1  => 8.31549
``````

Thank you for this.

1 Like

For future reference, here is the full example:

``````#Parameters
kb1 =  0.01;
K1  =  50.0;
d1  =  1.0;
l11 =  1.0;
m11 =  10.0;
l12 =  1.0;
m12 =  1.0;
kb2 =  0.01;
K2  =  20.0;
d2  =  1.0;
l22 =  5.0;
m22 =  1.0;
l21 =  10.0;
m21 =  1.0;

using JuMP
using Gurobi

model = JuMP.Model(Gurobi.Optimizer)
JuMP.set_optimizer_attribute(model, "NonConvex", 2)
JuMP.set_optimizer_attribute(model, "PoolSearchMode", 2)
JuMP.set_optimizer_attribute(model, "PoolSolutions", 100)

V = @variables(model, begin
X10
M1 >= 0
X11
M2 >= 0
X12
X20
X22
X21
end)

C = @constraints(model, begin
m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0
l11*M1*X10 - m11*X11 == 0
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0
l12*M2*X10 - m12*X12 == 0
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0
l22*M2*X20 - m22*X22 == 0
l21*M1*X20 - m21*X21 == 0
X10 + X12 + X11 == 1
X20 + X21 + X22 == 1
end)

optimize!(model)

#Seems to find only 1 fixed point M1=>40, M2=>0
result_dict = Dict(V .=> value.(V))

result_dicts = [Dict(V .=> value.(V; result=i)) for i in 1:result_count(model)]

"Determine if values of two `Dict`s compared on their common keys are all within `atol`"
function diff_dicts(d1::Dict, d2::Dict; atol=1e-4)

different_keys = []

for k in intersect(keys(d1), keys(d2))
if abs(d1[k] - d2[k]) > atol
push!(different_keys, k)
end
end

return isempty(different_keys)
end

candidates = [result_dicts[1]]

for dc1 in result_dicts
if dc1 in candidates
continue
end
if all([!diff_dicts(dc1, dc2) for dc2 in candidates])
push!(candidates, dc1)
end
end

display(candidates) # gives 3 different solutions
``````
2 Likes