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
- rewrite,
- or if you really insist, recover A x = b from a few evaluations of
F(egb = 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.
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:
[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};
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:
[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