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 [1] - x [2]
  f2 = 8 - x [1] - 3 * x [2]
  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] = 1 - x[1] - x[2]
         F[2] = 8 - x[1] - 3x[2]
       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[1] - x[2], 8 - x[1] - 3*x[2]]
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.

PS. Please quote your code.

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 :stuck_out_tongue_winking_eye:

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:
 [1] top-level scope at REPL[177]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088
G(x) = [3x[1] + 5x[1] + 2x[3] - 12, 2x[1] + 5x[2] - 4x[3] - 9, 4x[1] - x[2] + 5x[3] + 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
link

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[1] instead of 5x[2]

5 Likes

sys := { 3x + 5y + 2z = 12, 2x + 5y - 4z = 9,
4x - y + 5z = -3};

sys := {3x+5y+2z = 12, 2x+5y-4z = 9, 4x-y+5...

solve( sys, {x,y,z});

{y = 99/37, z = 24/37, x = -33/37}
Maple will automatically use fractions. However, you can force decimal answers using the evalf command.

evalf(%);

{y = 2.675675676, z = .6486486486, x = -.8918918919...

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:
 [1] top-level scope at REPL[189]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088

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