# Julia equivalent for python scipy.optimize.fsolve

What would be the Julia equivalent for python scipy.optimize.fsolve,
to resolve

``````function F (x)
f1 = 1- x  - x 
f2 = 8 - x  - 3 * x 
return f1, f2
end

x0 = [1,1]
``````

Solves a linear system of equations set in not matrix notation

1 Like

Have you tried a simple search for “optimization in Julia”? There are so many packages out there that it is hard to suggest just one solution. Start looking at Optim.jl, BlackBoxOptim.jl, …

3 Likes

What do you mean “set in not matrix notation”. If you have a function like `F` just pass `x -> collect(F(x))` to the solver, if the solver expects an array to be returned. Or do you mean it is a general non linear system?

2 Likes

If your system is linear, you really want `\`.

You can either

1. rewrite,
2. or if you really insist, recover A x = b from a few evaluations of `F` (eg `b = F(zeros(2))`, etc).
1 Like

There a several options, I think, but the NLsolve.jl package is one possibility:

``````julia> using NLsolve

julia> function F!(F, x)
F = 1 - x - x
F = 8 - x - 3x
end

julia> result = nlsolve(F!, [1.0,1.0], autodiff=:forward)
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [1.0, 1.0]
* Zero: [-2.5, 3.5]
* Inf-norm of residuals: 0.000000
* Iterations: 2
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 3
* Jacobian Calls (df/dx): 3

julia> result.zero
2-element Vector{Float64}:
-2.5
3.5
``````

Of course, if you know that your equations are linear, you should really put them into some kind of matrix and use `\`:

``````julia> [1 1; 1 3] \ [1, 8]
2-element Vector{Float64}:
-2.5
3.5
``````

but perhaps your real problem is nonlinear and you intended this `F` as a toy problem to figure out the library?

6 Likes

A. Solving Linear Systems

Here is a simple example which we can solve quite easily using the solve command.
sys := { 3x + 5y + 2z = 12, 2x + 5y - 4z = 9,4x - y + 5z = -3};
solve( sys, {x,y,z});

Or Exist in Julia a function or package that Convert a System of Linear Equations to Matrix Form?
Given a system of linear equations, determine the associated augmented matrix

Augmented Matrix for a Linear System
linear equations:

[3x + 5y + 2z = 12, 2x + 5y - 4z = 9,4x - y + 5z = -3]
List of variables:

[x,y,z]

Augmented matrix:
LinearAlgebraGenerateMatrix
3 5 2 12
2 5 -4 9
4 -1 5 -3

Yes, if you use `NLsolve` as in my post above with `autodiff=:forward` then it uses automatic differentiation to form the Jacobian matrix for Newton iterations. If you have a linear system of equations, the Jacobian matrix is exactly the matrix form of the problem, and Newton iterations will converge in a single `\` step.

And if you know that your equations are linear, you can do the differentiation directly using the ForwardDiff.jl package:

``````julia> using ForwardDiff

julia> F(x) = [1- x - x, 8 - x - 3*x]
F (generic function with 1 method)

julia> J = ForwardDiff.jacobian(F, [0,0])
2×2 Array{Int64,2}:
-1  -1
-1  -3

julia> x = -J \ F([0,0])
2-element Array{Float64,1}:
-2.5
3.5
``````

Here, the Jacobian `J` corresponds to the “matrix form” of your problem for a right-hand-side of `-F([0,0])`, and appropriate use of `\` (equivalent to a single Newton step) gives the same solution `x` as above.

7 Likes

ModellingToolkit.jl also works quite well for this kind of manipulation:

``````julia> using ModelingToolkit

julia> @variables x y z
(x, y, z)

julia> eqs = [3x + 5y + 2z - 12, 2x + 5y - 4z - 9, 4x - y + 5z + 3]
3-element Array{Operation,1}:
((3x + 5y) + 2z) - 12
((2x + 5y) - 4z) - 9
((4x - y) + 5z) + 3

julia> M = get.(ModelingToolkit.jacobian(eqs, [x, y, z]))
3×3 Array{Int64,2}:
3   5   2
2   5  -4
4  -1   5

julia> b = get.(substitute.(eqs, Ref([x, y, z] .=> 0)))
3-element Array{Int64,1}:
-12
-9
3

julia> M \ -b
3-element Array{Float64,1}:
-0.8918918918918918
2.675675675675676
0.6486486486486486

julia> substitute.(eqs, Ref([x, y, z] .=> ans))
3-element Array{ModelingToolkit.Constant,1}:
ModelingToolkit.Constant(0.0)
ModelingToolkit.Constant(0.0)
ModelingToolkit.Constant(0.0)
``````

To explicitely get the augmented matrix, you can just do the following:

``````julia> [M -b]
3×4 Array{Int64,2}:
3   5   2  12
2   5  -4   9
4  -1   5  -3
``````
5 Likes

Seeing lots of great options. I didn’t JuMP mentioned yet for solving larger linear programs, so here it is:

``````using JuMP, GLPK

m = model(GLPK.Optimizer)
@variable(m, x)
@variable(m, y)
@variable(m, z)
@constraints(m, begin
3x + 5y + 2z == 12,
2x + 5y - 4z == 9,
4x - y + 5z == -3
end)
optimize!(m)
value.((x,y,z))
``````

Where there above is just about the most verbose way to write it, but oh well 2 Likes
``````Julia>using JuMP, GLPK
[ Info: Precompiling JuMP [4076af6c-e467-56ae-b986-b466b2749572]
WARNING: using JuMP.@variables in module Main conflicts with an existing identifier.
[ Info: Precompiling GLPK [60bf3e95-4087-53dc-ae20-288a0d20c6a6]

Julia>m = model(GLPK.Optimizer)
ERROR: UndefVarError: model not defined
Stacktrace:
 top-level scope at REPL:1
``````
``````G(x) = [3x + 5x + 2x - 12, 2x + 5x - 4x - 9, 4x - x + 5x + 3]
J = ForwardDiff.jacobian(G, [0,0,0])
x = -J \ G([0,0,0])
7:36:28->>3-element Array{Float64,1}:
2.129032258064516
-1.064516129032258
-2.5161290322580645
``````

The expected results are:
x=-33/37 => -0.8918918918918919
y=99/37 => 2.675675675675676
z=24/37 =>0.6486486486486487

Why? `G(x)` gives you zero and using `x = [-33/37, 99/37, 24/37]` does not give a zero.

apologies. It should be `Model`, not `model`

In MATLAB, this is done in a single line of code, I also wonder if there is a similar thing in Julia:

``````>> syms x y z
>> [x,y,z] = solve(3*x + 5*y + 2*z == 12, 2*x + 5*y - 4*z == 9,4*x - y + 5*z == -3)

x =
-33/37

y =
99/37

z =
24/37
``````

the result shown in

``````M \ -b
3-element Array{Float64,1}:
-0.8918918918918918
2.675675675675676
0.6486486486486486
``````

The expected results are:
x=-33/37 => -0.8918918918918919
y=99/37 => 2.675675675675676
z=24/37 =>0.6486486486486487

You have `5x` instead of `5x`

5 Likes

sys := { 3x + 5y + 2z = 12, 2x + 5y - 4z = 9,
4x - y + 5z = -3}; solve( sys, {x,y,z}); Maple will automatically use fractions. However, you can force decimal answers using the evalf command.

evalf(%); ``````Julia>m = Model(GLPK.Optimizer)
ERROR: MethodError: no method matching Model(::Type{GLPK.Optimizer})
Closest candidates are:
Model(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) at C:\Users\hermesr\.julia\packages\JuMP\iGamg\src\JuMP.jl:126
Model(; caching_mode, solver) at C:\Users\hermesr\.julia\packages\JuMP\iGamg\src\JuMP.jl:161
Model(::MathOptInterface.AbstractOptimizer, ::Dict{MathOptInterface.ConstraintIndex,AbstractShape}, ::Set{Any}, ::Any, ::Any, ::Dict{Symbol,Any}, ::Int64, ::Dict{Symbol,Any}) at C:\Users\hermesr\.julia\packages\JuMP\iGamg\src\JuMP.jl:126
...
Stacktrace:
 top-level scope at REPL:1
``````

Yes, in MATLAB, you use `double([x y z])` to get decimal answer.

There are various packages for symbolic manipulations, but that’s probably overkill for this application.

I am curious about the context here. If one is talking about a single linear system, it is trivial to transcribe it to a matrix form. If, on the other hand, one needs to solve a lot of linear systems, it is hard to imagine that they come in the general form and not as a matrix equation.

``````julia> A = [3 5 2;
2 5 -4;
4 -1 5];

julia> b = [12, 9, -3];

julia> Rational.(A) \ Rational.(b)
3-element Array{Rational{Int64},1}:
-33//37
99//37
24//37
``````

You may want to use `Rational{BigInt}` in general though.

That said, I am wondering if you misunderstand what Julia is for: while various packages address this, Julia in general is not a CAS like Maple, Mathematica, or Maxima.

5 Likes