# Numeric value substitution into SymPy Symbolic array or matrix

**URL:** <https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204>\
**Category:** Numerics\
**Tags:** sympy\
**Created:** [May 22, 2023, 8:34am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204 "2023-05-22T08:34:51Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [May 22, 2023, 8:34am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204/1 "2023-05-22T08:34:51Z")

</div>

I am having trouble substituting numerical values into symbolic arrays and matrices.  
I have read about [`lambdify`](https://docs.juliahub.com/SymPy/KzewI/1.1.8/Tutorial/basic_operations/#lambdify-1), but I still fail to provide numerical values for the arguments in a consistent way.

As an example, consider the following `Julia` code

```julia
julia> using Sympy
julia> n = 4;
julia> η = [symbols("η$i", positive=true) for i in 1:n];
julia> α = [(i<3) ? 1 : symbols("α", positive=true) for i in 1:n];
julia> η .*= α
4-element Vector{Sym}:
   η₁
   η₂
 α⋅η₃
 α⋅η₄

julia> syms = Symbol.(free_symbols(tensors[:η]))
5-element Vector{Symbol}:
 :α
 :η1
 :η2
 :η3
 :η4

julia> f = lambdify(η, syms)
#122 (generic function with 1 method)

julia> f.(0.5,1.0,1.0,1.0,1.0)
4-element Vector{Float64}:
 1.0
 1.0
 0.5
 0.5

```

As expected, the outcome is indeed the substitution for `α=0.5` and all `η=1.0`. However, my problem arises from the fact that I would like to feed the expressions from arrays or tuples, i.e. what I would like is:

```julia
julia> values = (0.5, 1.0, 1.0, 1.0, 1.0);
julia> f.(values)
ERROR: MethodError: no method matching var"##339"(::Float64)

Closest candidates are:
  var"##339"(::Any, ::Any, ::Any, ::Any, ::Any)
   @ SymPy none:0

Stacktrace:
  [1] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:816
  [2] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base ./essentials.jl:813
  [3] (::SymPy.var"#118#119"{SymPy.var"###339"})(::Float64, ::Vararg{Float64})
    @ SymPy ~/.julia/packages/SymPy/gkr7m/src/lambdify.jl:223
  [4] map(::Function, ::Float64)
    @ Base ./number.jl:284
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
  [7] getindex
    @ ./broadcast.jl:610 [inlined]
  [8] copy
    @ ./broadcast.jl:912 [inlined]
  [9] materialize
    @ ./broadcast.jl:873 [inlined]
 [10] (::SymPy.var"#122#124"{Vector{SymPy.var"#118#119"}})(args::Float64)
    @ SymPy ~/.julia/packages/SymPy/gkr7m/src/lambdify.jl:253
 [11] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
 [12] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
 [13] #31
    @ ./broadcast.jl:1088 [inlined]
 [14] macro expansion
    @ ./ntuple.jl:72 [inlined]
 [15] ntuple
    @ ./ntuple.jl:69 [inlined]
 [16] copy
    @ ./broadcast.jl:1088 [inlined]
 [17] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple}, Nothing, SymPy.var"#122#124"{Vector{SymPy.var"#118#119"}}, Tuple{NTuple{5, Float64}}})
    @ Base.Broadcast ./broadcast.jl:873
 [18] top-level scope
    @ REPL[199]:1

```

The reason that I want this is so that I can make all of my symbolic arrays into “lambdified” expressions and evaluate them, for example at multiple points in parallel. Additionally, when the array is instead a matrix, I also do not wish to manually input all the matrix elements – especially when the sizes of the matrices grow.  
Consider for example,

```julia
julia> M = [symbols("M$i$j") for i in 1:3 for j in 1:3].reshape(3,3)
3×3 Matrix{Sym}:
 M₁₁ M₁₂ M₁₃
 M₂₁ M₂₂ M₂₃
 M₃₁ M₃₂ M₃₃
julia> syms_M = free_symbols.(M)
3×3 Matrix{Vector{Sym}}:
 [M11] [M12] [M13]
 [M21] [M22] [M23]
 [M31] [M32] [M33]
julia> fM = lambdify(M, syms_M)
#122 (generic function with 1 method)

```

How do I use `fM` to return a numeric valued matrix, e.g. when `Mvals = rand(3,3)`? Any help would be greatly appreciated!

---

<div class="post-metadata">

**Author:** ![sijo](https://avatars.discourse-cdn.com/v4/letter/s/da6949/32.png) [@sijo](https://discourse.julialang.org/u/sijo)\
**Post date:** [May 22, 2023, 9:12am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204/2 "2023-05-22T09:12:27Z")

</div>

You would write

```julia
julia> values = (0.5, 1.0, 1.0, 1.0, 1.0);
julia> f.(values)

```

if `f` is a function of one parameter and you want to call it 5 times and gather the 5 results in a vector.

In your case however `f` is a function that takes 5 parameters. So to call it with a tuple of values you can write

```julia
julia> f(values...)
4-element Vector{Float64}:
 1.0
 1.0
 0.5
 0.5

```

For the matrix example: here you broadcast the `free_symbols` call. I don’t think that will work: it returns a matrix where each element is the result of a `free_symbols` call, i.e. each element is a vector of symbols. When you call `lambdify` with this matrix of vectors of symbols, for some reason it returns something instead of raising an error, but the returned function doesn’t seem usable. It’s not surprising since you should call `lambdify` with a vector of symbols instead.

To solve this you could convert the matrix of vectors into a simple vector. But the easiest is to call `free_symbols` directly with a vector:

```julia
julia> M = [symbols("M$i$j") for i in 1:3, j in 1:3]
3×3 Matrix{Sym}:
 M₁₁ M₁₂ M₁₃
 M₂₁ M₂₂ M₂₃
 M₃₁ M₃₂ M₃₃

julia> syms_M = free_symbols(vec(M))
9-element Vector{Sym}:
 M₁₁
 M₁₂
 M₁₃
 M₂₁
 M₂₂
 M₂₃
 M₃₁
 M₃₂
 M₃₃

julia> fM = lambdify(M, syms_M)
#122 (generic function with 1 method)

julia> Mvals = rand(3,3)
3×3 Matrix{Float64}:
 0.0253515 0.261792 0.640684
 0.668682 0.216239 0.433714
 0.780887 0.808188 0.315145

julia> fM(Mvals...)
3×3 Matrix{Float64}:
 0.0253515 0.668682 0.780887
 0.261792 0.216239 0.808188
 0.640684 0.433714 0.315145

```

Note: I think the splatting of a matrix with `...` is not very efficient. If that’s a problem, maybe you can solve it by using a matrix type from StaticArrays.jl (so that the matrix size is known by the compiler).

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [May 22, 2023, 9:30am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204/3 "2023-05-22T09:30:52Z")

</div>

Thanks for the quick response. I was just starting to type up this, as I just figured it out. It makes indeed sense to use splatting. The efficiency of the splatting is not really a deal breaker, but I might need to convert to a `StaticArray` if it appears to be a problem.

Another question I have is that I do not know how to assign values to elements that are a function of multiple variables. I can do it in this way, most likely, but then I need to know the order of symbols that `free_symbols(expr)` returns – which is odd as I do not necessarily know this beforehand.

Consider for example a simple Lotka-Volterra system (which are similar to the systems I need this substitution for):

```julia
julia> M = reshape([-symbols("a$i$j", positive=true) for i in 1:2 for j in 1:2], (2,2))
2×2 Matrix{Sym}:
 -a₁₁ -a₂₁
 -a₁₂ -a₂₂

julia> r = diagm([symbols("η", positive=true), -symbols("μ", positive=true)])
2×2 Matrix{Sym}:
 η 0
 0 -μ

julia> A = r+M
2×2 Matrix{Sym}:
 -a₁₁ + η -a₂₁
     -a₁₂ -a₂₂ - μ

julia> syms = free_symbols(vec(A))
4-element Vector{Sym}:
 a₁₂
 a₂₁
   η
   μ

julia> Avals = vec([0 1.0; 1.0 0]);
julia> ηval = 1.0;
julia> μval = 0.5;
julia> values = vcat(Avals, ηval, μval);
julia> f = lambdify(A, syms)
julia> f(values...)
2×2 Matrix{Float64}:
  1.0 -1.0
 -1.0 -0.5

```

as required. However, this forces me to know the order of variables that `free_symbols(expr)` returns, as if I instead have `values = vcat(ηval, μval, Avals)`, then `f(values...)` will understandably be wrong. Is there a way to have `lambdify()` return a function with named arguments, or something like that? Or is there another way that I can safely insert numerical values of each symbol subsequently, so that I know for sure that substituting the value of a symbol results in the evaluation of only that symbol within the expression? I hope my question is clear.

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [May 22, 2023, 10:36am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204/4 "2023-05-22T10:36:22Z")

</div>

So diving into the source code of `SymPy` we can find the [`free_symbols()` function definition](https://docs.juliahub.com/SymPy/KzewI/1.1.8/reference/#SymPy.free_symbols-Union%7BTuple%7BUnion%7BVector%7BT%7D,%20T%7D%7D,%20Tuple%7BT%7D%7D%20where%20T%3C:SymPy.SymbolicObject). In there, one can see that the sorting is executed by calling `sortperm`:

```julia
fs[sortperm(string.(fs))] # some sorting order to rely on

```

This gives us a way to reliably generate the order, without having to call it beforehand. Using this we are free to specify a dictionary that holds all the values corresponding to the symbols present in the (matrix) expression, order the accordingly, and call the `lambdify`’d function with the appropriate splatted argument list.  
For example:

```julia
julia> using SymPy, LinearAlgebra

julia> A = reshape([-symbols("a$i$j", positive=true) for i in 1:2 for j in 1:2], (2,2))
2×2 Matrix{Sym}:
 -a₁₁ -a₂₁
 -a₁₂ -a₂₂

julia> r = diagm([symbols("η", positive=true), -symbols("μ", positive=true)])
2×2 Matrix{Sym}:
 η 0
 0 -μ

julia> Ω = r + A
2×2 Matrix{Sym}:
 -a₁₁ + η -a₂₁
     -a₁₂ -a₂₂ - μ

julia> Avals = [0.0 1.0; 1.0 0.0]
2×2 Matrix{Float64}:
 0.0 1.0
 1.0 0.0

julia> adict = Dict([Pair(symbols("a$i$j"), Avals[i,j]) for i in 1:2 for j in 1:2])
Dict{Sym, Float64} with 4 entries:
  a22 => 0.0
  a11 => 0.0
  a12 => 1.0
  a21 => 1.0

julia> ηdict = Dict(Pair(symbols("η"), 1.0))
Dict{Sym, Float64} with 1 entry:
  η => 1.0

julia> μdict = Dict(Pair(symbols("μ"), 0.5))
Dict{Sym, Float64} with 1 entry:
  μ => 0.5

julia> var2val = merge(adict, ηdict, μdict)
Dict{Sym, Float64} with 6 entries:
  a22 => 0.0
  η => 1.0
  μ => 0.5
  a11 => 0.0
  a12 => 1.0
  a21 => 1.0

julia> order = sortperm(string.(keys(var2val)))
6-element Vector{Int64}:
 4
 5
 6
 1
 2
 3

julia> vals = collect(values(var2val))[order]
6-element Vector{Float64}:
 0.0
 1.0
 1.0
 0.0
 1.0
 0.5

julia> fΩ = lambdify(Ω, free_symbols(vec(Ω)))
#122 (generic function with 1 method)

julia> fΩ(vals...)
2×2 Matrix{Float64}:
  1.0 -1.0
 -1.0 -0.5

```

This can probably be done more cleanly considering the available parameters, but the one that defines the problem will probably know the parameters and can therefore construct the `var2val` dictionary without any issue.

---

<div class="post-metadata">

**Author:** ![sijo](https://avatars.discourse-cdn.com/v4/letter/s/da6949/32.png) [@sijo](https://discourse.julialang.org/u/sijo)\
**Post date:** [May 22, 2023, 10:39am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204/5 "2023-05-22T10:39:50Z")

</div>

The `free_symbols` function is just a convenient way to get a list of symbols, but you don’t _have to_ use it. If you know the symbols you can build the list yourself explicitly, and this way you control the order:

```julia
julia> A0 = [symbols("a$i$j", positive=true) for i in 1:2, j in 1:2]
2×2 Matrix{Sym}:
 a₁₁ a₁₂
 a₂₁ a₂₂

julia> @vars η μ positive=true
(η, μ)

julia> r = [η 0; 0 μ];

julia> A = r-A0
2×2 Matrix{Sym}:
 -a₁₁ + η -a₁₂
     -a₂₁ -a₂₂ + μ

julia> Avals = vec([0 1.0; 1.0 0]);

julia> ηval = 1.0;

julia> μval = 0.5;

julia> syms = [vec(A0); η; μ] # here I choose an order for the symbols
6-element Vector{Sym}:
 a₁₁
 a₂₁
 a₁₂
 a₂₂
   η
   μ

julia> f = lambdify(A, syms);

julia> f(values...)
2×2 Matrix{Float64}:
  1.0 -1.0
 -1.0 0.5

```

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [May 22, 2023, 10:42am UTC](https://discourse.julialang.org/t/numeric-value-substitution-into-sympy-symbolic-array-or-matrix/99204/6 "2023-05-22T10:42:10Z")

</div>

That makes sense as well… and is probably a bit easier. Thanks!