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

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`

(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 := { 3*x + 5*y + 2*z = 12, 2*x + 5*y - 4*z = 9,4*x - y + 5*z = -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:

[3*x + 5*y + 2*z = 12, 2*x + 5*y - 4*z = 9,4*x - y + 5*z = -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 := { 3*x + 5*y + 2*z = 12, 2*x + 5*y - 4*z = 9,

4*x - y + 5*z = -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