Day_07 - Julia’s Symbolics: Computer Algebra System

Day_07 - Julia’s `Symbolics`: Computer Algebra System#

Let the computer do the hard work#

The double-pendulum from Day_06 was a great success. I got some experience building differential equations, pseudo-benchmarking some results, and building some excellent animations. One of the tough parts of physical simulations is getting all the math worked out correctly. Its satisfying to eliminate terms and come up with a final equation of motion, kind of like a Sudoku puzzle or Wordle game, but it takes time. When you start the next problem, you’re back to square one.

Enter Computer Algebra Systems #

Computer Algebra Systems (CAS) allow you to create objects that follow symbolic algebra, calculus, etc. I avoided CAS programs for years because it felt so difficult to get the equations back into a numerical form. As an engineer, I mostly deal with tables, fractions, percents, etc. I don’t always need or want an exact function. That’s why CFD, FEA, and multibody dynamics are so useful to engineers.

Julia’s CAS `using Symbolics`#

There are some great examples on Symbolics.jl documentation. Here, I want to

build some simple sets of functions
take some derivatives
create a callable numerical function (not a CAS object)

using Plots, DifferentialEquations, LaTeXStrings, Symbolics

[ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] (cache misses: wrong dep version loaded (2))

Failed to precompile Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80] to "/home/ryan/.julia/compiled/v1.11/Plots/jl_HynGIh".

Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool; flags::Cmd, cacheflags::Base.CacheFlags, reasons::Dict{String, Int64}, isext::Bool)
    @ Base ./loading.jl:3085
  [3] (::Base.var"#1082#1083"{Base.PkgId})()
    @ Base ./loading.jl:2492
  [4] mkpidlock(f::Base.var"#1082#1083"{Base.PkgId}, at::String, pid::Int32; kwopts::@Kwargs{stale_age::Int64, wait::Bool})
    @ FileWatching.Pidfile ~/projects/julia/usr/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:95
  [5] #mkpidlock#6
    @ ~/projects/julia/usr/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:90 [inlined]
  [6] trymkpidlock(::Function, ::Vararg{Any}; kwargs::@Kwargs{stale_age::Int64})
    @ FileWatching.Pidfile ~/projects/julia/usr/share/julia/stdlib/v1.11/FileWatching/src/pidfile.jl:116
  [7] #invokelatest#2
    @ ./essentials.jl:1057 [inlined]
  [8] invokelatest
    @ ./essentials.jl:1052 [inlined]
  [9] maybe_cachefile_lock(f::Base.var"#1082#1083"{Base.PkgId}, pkg::Base.PkgId, srcpath::String; stale_age::Int64)
    @ Base ./loading.jl:3609
 [10] maybe_cachefile_lock
    @ ./loading.jl:3606 [inlined]
 [11] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:2488
 [12] __require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2315
 [13] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
 [14] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
 [15] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2302
 [16] macro expansion
    @ ./loading.jl:2241 [inlined]
 [17] macro expansion
    @ ./lock.jl:273 [inlined]
 [18] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2198
 [19] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
 [20] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
 [21] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2191

Build a set of functions and the Jacobian#

The first thing I build is a set of functions and its Jacobian,

$\mathbf{f}(x_1,x_2) = \left[ \begin{array}{c} x_2 + x_1^{2} \ x_1 + x_2^{2} \ \end{array} \right] $

$\frac{\partial f_i}{\partial x_j} = \left[ \begin{array}{cc} 2 x_{1} & 1 \ 1 & 2 x_{2} \ \end{array} \right]$

There were a few tricks to get everything in the right place and cooperating. So here is the process

define the variables x_1 and x_2
define $\mathbf{f(x)}$ f
build the function f_x with build_function
evaluate the first returned output from build_function

1. define the variables `x_1` and `x_2`#

@variables x_1, x_2

\[\begin{split}\begin{equation} \left[ \begin{array}{c} x_{1} \\ x_{2} \\ \end{array} \right] \end{equation}\end{split}\]

2. define $\mathbf{f(x)}$ `f`#

Here, I use the same array definitions as Julia to build a vector of functions.

f = [x_1^2 + x_2, x_2^2 + x_1]

\[\begin{split}\begin{equation} \left[ \begin{array}{c} x_{2} + x_{1}^{2} \\ x_{1} + x_{2}^{2} \\ \end{array} \right] \end{equation}\end{split}\]

3. build the function `f_x` with `build_function`#

To get a reusable, numerical function, I want to define a Julia function. I could hard code it in, but then I’m back to doing all the math myself. The build_function is a fast Julia interpretation of the algebra contained in the CAS function. It returns a full function that can be saved as a *.jl file.

Note: There is a new term to me called the world-age issues that has to do with where and when the function is defined. This method using eval creates a new world within the Julia session that cannot be precompiled. I believe saving this to the *.jl file is a workaround to this slow-down.

f_x= build_function(f, x_1, x_2, expression = Val{false})

(RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:x_1, :x_2), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x6a792526, 0xcbf59f6b, 0x6d67a95f, 0x96e674fe, 0x64d4827f)}(quote
    #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:349 =#
    #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:350 =#
    #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:351 =#
    begin
        #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:444 =#
        (SymbolicUtils.Code.create_array)(Array, nothing, Val{1}(), Val{(2,)}(), (+)(x_2, (^)(x_1, 2)), (+)(x_1, (^)(x_2, 2)))
    end
end), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :x_1, :x_2), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x184abcc3, 0x978c6f1e, 0x2ab512ab, 0xe6b23c0f, 0xaf4ffa04)}(quote
    #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:349 =#
    #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:350 =#
    #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:351 =#
    begin
        #= /home/ryan/.julia/packages/Symbolics/1OrKJ/src/build_function.jl:378 =#
        #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:398 =# @inbounds begin
                #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:394 =#
                ˍ₋out[1] = (+)(x_2, (^)(x_1, 2))
                ˍ₋out[2] = (+)(x_1, (^)(x_2, 2))
                #= /home/ryan/.julia/packages/SymbolicUtils/v2ZkM/src/code.jl:396 =#
                nothing
            end
    end
end))

4. evaluate the first returned output from `build_function`#

In this case, there are two functions:

The first is a function f([x, y]) that computes and builds an output vector $[x_1^2 + x_2,~x_2^2 + x_1]$.
The second function is an in-place non-allocating mutating function which mutates its first value

I don’t need the second function at the moment, so I evaluate the first function.

out = f_x[1](1.0, 2.0)

2-element Vector{Float64}:
 3.0
 5.0

Plots.theme(:cooper)
x1_array = range(0, 1, length = 100)
x2_array = range(0, 1, length = 100)
Farray1 = [f_x[1](x1i, x2i)[1] for x1i in x1_array for x2i in x2_array]
Farray2 = [f_x[1](x1i, x2i)[2] for x1i in x1_array for x2i in x2_array]

l = @layout [a; b]
p = contourf(x1_array, x2_array, Farray1, layout = l)
contourf!(p[2], x1_array, x2_array, Farray2)

_images/755c075e35f05f10e3383a082d89cbc7ba2c42e15ddf8c76928f7c1d69638b75.svg

Evaluate the Jacobian#

Success! Now, I can calculate the Jacobian $\frac{\partial f_i}{\partial x_i}$ and evaluate its result with the same steps.

First, create the symbolic function, DJ, with .jacobian.

DJ = Symbolics.jacobian(f, [x_1, x_2])

\[\begin{split}\begin{equation} \left[ \begin{array}{cc} 2 x_{1} & 1 \\ 1 & 2 x_{2} \\ \end{array} \right] \end{equation}\end{split}\]

Next, build_function and eval the result.

DJ_xy = build_function(DJ, [x_1, x_2], expression = Val{false})
DJ_xy[1]([3, 2])

2×2 Matrix{Int64}:
 6  1
 1  4

Success!

Make a time-varying array of functions#

Now, I want to create a set of functions that vary with time, take the derivative and evaluate the results.

This will give me the piece I need for the $\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right)$ calculations.

I want to build,

$\mathbf{g(z)} = \left[ \begin{array}{c} \mathrm{z{_1}}\left( t \right) \mathrm{z{_2}}\left( t \right) \ 2 \mathrm{z{_4}}\left( t \right) + \mathrm{z{_2}}\left( t \right) + \mathrm{z{_3}}\left( t \right) \

\left( \mathrm{z{_4}}\left( t \right) \right)^{2} + \left( \mathrm{z{_1}}\left( t \right) \right)^{2} \ \end{array} \right]$

then, take the derivative,

$\frac{d}{dt}(\mathbf{g(z)})=\left[ \begin{array}{c} \mathrm{z{_1}}\left( t \right) \frac{dz{_2}(t)}{dt} + \mathrm{z{_2}}\left( t \right) \frac{dz{_1}(t)}{dt} \ 2 \frac{dz{_4}(t)}{dt} + \frac{dz{_2}(t)}{dt} + \frac{dz{_3}(t)}{dt} \ 2 \mathrm{z{_1}}\left( t \right) \frac{dz{_1}(t)}{dt} - 2 \mathrm{z{_4}}\left( t \right) \frac{dz{_4}(t)}{dt} \ \end{array} \right]$

Build a $d/dt$ operator#

The next function, I make is just an operator. It takes derivatives of variables that are a function of t, time. Now, I can create differential equations.

First, define t

@variables t

\[\begin{split}\begin{equation} \left[ \begin{array}{c} t \\ \end{array} \right] \end{equation}\end{split}\]

Now, define the $\frac{d}{dt}$ as Differential(t)

D = Differential(t)

(::Differential) (generic function with 2 methods)

Now, I create four variables, $\mathbf{z} = [z_1,~z_2,~z_3,z_4]$ that are functions of time. I use the shortcut to create an array of symbols z[1:4] and the (t) to indicate that each variable is a function of time.

@variables z[1:4](t)
z = [zi for zi in z]

\[\begin{split}\begin{equation} \left[ \begin{array}{c} \mathrm{z{_1}}\left( t \right) \\ \mathrm{z{_2}}\left( t \right) \\ \mathrm{z{_3}}\left( t \right) \\ \mathrm{z{_4}}\left( t \right) \\ \end{array} \right] \end{equation}\end{split}\]

Then, I can define the derivatives of each of these variables.

dz = D.(z)

\[\begin{split}\begin{equation} \left[ \begin{array}{c} \frac{dz{_1}(t)}{dt} \\ \frac{dz{_2}(t)}{dt} \\ \frac{dz{_3}(t)}{dt} \\ \frac{dz{_4}(t)}{dt} \\ \end{array} \right] \end{equation}\end{split}\]

Next, build the function $\mathbf{g(z)}$

g = [z[1]*z[2], z[3]+z[2]+2*z[4], z[1]^2-z[4]^2]

\[\begin{split}\begin{equation} \left[ \begin{array}{c} \mathrm{z{_1}}\left( t \right) \mathrm{z{_2}}\left( t \right) \\ 2 \mathrm{z{_4}}\left( t \right) + \mathrm{z{_2}}\left( t \right) + \mathrm{z{_3}}\left( t \right) \\ - \left( \mathrm{z{_4}}\left( t \right) \right)^{2} + \left( \mathrm{z{_1}}\left( t \right) \right)^{2} \\ \end{array} \right] \end{equation}\end{split}\]

When I take the derivative of g with D., it is a lazy evaluation. To get the expanded result, I call out expand_derivatives from Symbolics.jl.

D.(g)

\[\begin{split}\begin{equation} \left[ \begin{array}{c} \mathrm{\frac{d}{d t}}\left( \mathrm{z{_1}}\left( t \right) \mathrm{z{_2}}\left( t \right) \right) \\ \mathrm{\frac{d}{d t}}\left( 2 \mathrm{z{_4}}\left( t \right) + \mathrm{z{_2}}\left( t \right) + \mathrm{z{_3}}\left( t \right) \right) \\ \mathrm{\frac{d}{d t}}\left( - \left( \mathrm{z{_4}}\left( t \right) \right)^{2} + \left( \mathrm{z{_1}}\left( t \right) \right)^{2} \right) \\ \end{array} \right] \end{equation}\end{split}\]

Note: “Lazy” evaluation is not a negative comment on the function. I’m just pointing out that it doesn’t fully evaluate the result until I think it’s necessary. Its a feature not a bug.

Dg = expand_derivatives.(D.(g))

\[\begin{split}\begin{equation} \left[ \begin{array}{c} \mathrm{z{_1}}\left( t \right) \frac{dz{_2}(t)}{dt} + \mathrm{z{_2}}\left( t \right) \frac{dz{_1}(t)}{dt} \\ 2 \frac{dz{_4}(t)}{dt} + \frac{dz{_2}(t)}{dt} + \frac{dz{_3}(t)}{dt} \\ 2 \mathrm{z{_1}}\left( t \right) \frac{dz{_1}(t)}{dt} - 2 \mathrm{z{_4}}\left( t \right) \frac{dz{_4}(t)}{dt} \\ \end{array} \right] \end{equation}\end{split}\]

Tripping point for z(t) as a callable variable#

This was a wall that I hit, hopefully you can learn from my mistakes. When I try to build_function and eval the resulting code, I get this error that (z[1](t)) is not callable. I believe it has to do with the assumption that z is a function of time on my part.

julia> build_function(Dg, z, dz)
julia> eval(Dg[1])([1,3,5,7], [2,4,6,8])

Sym z[1](t)*Differential(t)(z[2](t)) + z[2](t)*Differential(t)(z[1](t)) is not callable. Use @syms z[1](t)*Differential(t)(z[2](t)) + z[2](t)*Differential(t)(z[1](t))(var1, var2,...) to create it as a callable.

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] (::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})(::Vector{Int64}, ::Vararg{Vector{Int64}})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/v2ZkM/src/types.jl:172
 [3] (::Num)(::Vector{Int64}, ::Vararg{Vector{Int64}})
   @ Symbolics ~/.julia/packages/Symbolics/HDE84/src/num.jl:17
 [4] top-level scope
   @ In[267]:2
 [5] eval
   @ ./boot.jl:373 [inlined]
 [6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196

Work-around for non-callable function variables#

Here is my workaround:

create callable variables q and dq with the same length as z(t)
substitute the values into the Dg symbolic function and simplify the results
build the function with the q and dq callable variables

I substitute the dq variables first, so I don’t end up with $\frac{d}{dt}(q_i)$-terms. The susbstitute function acts as such

substitute. makes substitutions on each function within the vector Dg
Dict(dz[i] => dq[i]) and Dict(z[i] => q[i]) replaces each $z_i(t)$ with $q_i$
for i in range(1,4) loops the dictionary from z[1] to z[4]

@variables q[1:4]
@variables dq[1:4]
Dgdq = simplify.(substitute.(Dg, (Dict(dz[i] => dq[i] for i in range(1, 4)),)))
gq = simplify.(substitute.(g, (Dict(z[i] => q[i] for i in range(1, 4)),)))
Dgq = simplify.(substitute.(Dgdq, (Dict(z[i] => q[i] for i in range(1, 4)),)))

\[\begin{split}\begin{equation} \left[ \begin{array}{c} dq{_2} q{_1} + dq{_1} q{_2} \\ 2 dq{_4} + dq{_2} + dq{_3} \\ 2 dq{_1} q{_1} - 2 dq{_4} q{_4} \\ \end{array} \right] \end{equation}\end{split}\]

Finally, I can build_function from Dgq and eval the result.

g_func = build_function(gq, q, expression = Val{false})
Dg_func = build_function(Dgq, q, dq, expression = Val{false})
Dg_func[1]([1,3,5,7], [2,4,6,8])
# eval(Dg_func[1])([1,3,5,7], [2,4,6,8])

3-element Vector{Int64}:
   10
   26
 -108

Animate some function evaluations#

Finally, I want to build an animation of the functions I am building, $\mathbf{g(z)}~and~\frac{d\mathbf{g}}{dt}$. I have to define my numerical variable space

time, t_space from 0 - $2\pi$ seconds
$\mathbf{z}$, z_space as $\sin(t+\phi)$
$\mathbf{\dot{z}}$, dz_space as $\cos(t+\phi)$
$\mathbf{g(z)}$, g_val as the evaluated g_func
$\mathbf{\dot{g}(z,\dot{z})}$, Dg_val as the evaluated Dg_func

t_space = range(0,2*pi, length = 200)

0.0:0.03157379551346526:6.283185307179586

Here, I plot the values for $z~and~\dot{z}$.

l = @layout [a; b]
z_space = zeros(length(t_space), 4)
dz_space = zeros(length(t_space), 4)
p = plot(title = "", layout = l)
for (i, phi) in enumerate([0, 0.2, 0.4, 0.6])
    z_space[:, i] = sin.(t_space .+ phi)
    dz_space[:, i] = cos.(t_space .+ phi)
    plot!(p[1], t_space, sin.(t_space .+ phi), label = "z_$i")
    plot!(p[2], t_space, cos.(t_space .+ phi), label = "dz_$i/dt")
    plot
end
p

_images/438e2671b0ba5e04a643693f614abd33d8d6c3914a5a72e9e9b7da3367040fb8.svg

Now, I plot the evaluated $\mathbf{g(z)}~and~\mathbf{\dot{g}(z)}$

l = @layout [a; b]
p = plot(xlabel = "time", layout = l, xlim = (0, 2*pi), ylim = (-5, 5))
g_val = zeros(length(t_space), 3)
Dg_val = zeros(length(t_space), 3)
for i in range(1, length(t_space))
    g_val[i, :] = g_func[1](z_space[i, :])
    Dg_val[i, :] = Dg_func[1](z_space[i, :], dz_space[i, :])
end
for i in range(1, 3)
    plot!(p[1], t_space, g_val[:, i], axis = nothing, label = "g_$i")
    plot!(p[2], t_space, Dg_val[:, i], axis = nothing, label = "dg_$i/dt")
end
p

_images/67fcd67fe7e8419208a3c44f9d763ac579bfaa40ca02fbadf170e6fbf2f3d192.svg

Finally, I can combine the input+evaluations into one animated gif to see the effects of $\mathbf{z}(t)$ on $\mathbf{\dot{z},~g,~\dot{g}}$.

l = @layout [a b; c d]
p = plot(xlabel = "time", layout = l, xlim = (0, 2*pi))#, ylim = (-5, 5))
import ColorSchemes.rainbow

@gif for ti in range(1, length(t_space))
    
    for i in range(1, 4)
        color = get(rainbow, (i-1)/(3))
        plot!(p[1], t_space[1:ti], z_space[1:ti, i], axis = nothing, label = "", title="z", linecolor = color)
        plot!(p[3], t_space[1:ti], dz_space[1:ti, i], axis = nothing, label = "", title = "dz/dt", linecolor = color)
    end
    for i in range(1, 3)
        color = get(rainbow, (i-1)/(2))
        plot!(p[2], t_space[1:ti], g_val[1:ti, i], axis = nothing, label = "", title = "g(z)",linecolor = color)
        plot!(p[4], t_space[1:ti], Dg_val[1:ti, i], axis = nothing, label = "", title = "dg/dt", linecolor = color)
    end
end

┌ Info: Saved animation to 
│   fn = /home/ryan/Documents/Career_docs/cooperrc-gh-pages/Julia-learning/tmp.gif
└ @ Plots /home/ryan/.julia/packages/Plots/D9pfj/src/animation.jl:114

print("Plotted above are the evaluations for g1(z1,z2,z3,z4), g2(z1,z2,z3,z4), g3(z1,z2,z3,z4)\n")
g

Plotted above are the evaluations for g1(z1,z2,z3,z4), g2(z1,z2,z3,z4), g3(z1,z2,z3,z4)

\[\begin{split}\begin{equation} \left[ \begin{array}{c} \mathrm{z{_1}}\left( t \right) \mathrm{z{_2}}\left( t \right) \\ 2 \mathrm{z{_4}}\left( t \right) + \mathrm{z{_2}}\left( t \right) + \mathrm{z{_3}}\left( t \right) \\ - \left( \mathrm{z{_4}}\left( t \right) \right)^{2} + \left( \mathrm{z{_1}}\left( t \right) \right)^{2} \\ \end{array} \right] \end{equation}\end{split}\]

Success!

Wrapping up#

Today I spent a chunk of time building and debugging symbolic algebra code. To me, the hardest part is always getting the CAS objects back into a function that can be evaluated, integrated, optimized, etc. The $\LaTeX$ output of the Symbolics.jl package are excellent. The ability to create differential operators and set up equations is amazing. I was able to define arrays of functions with multiple variables, take partial derivatives with the .jacobian, and full derivatives with Differential(...).

I am still wrapping my head around the output from build_function. It seems like the best practice is to create a seperate *.jl file, but for now I got the numerical functions to evaluate based upon their symbolic counterparts. The world-age problem is new to me. Coming from Python+MATLAB, I hadn’t really worried about or considered where and when my functions were compiled. I would guess that getting these concepts in Julia would also help in Python/MATLAB.

Coming up next: I’ll build a set of differential equations based upon the Lagrange equations $\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i}= 0$.

Day_07 - Julia’s Symbolics: Computer Algebra System

Contents

Day_07 - Julia’s Symbolics: Computer Algebra System#