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

Specific Domains
, , ,
4

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