# 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

1. build some simple sets of functions

2. take some derivatives

3. create a callable numerical function (not a CAS object)

using Plots, DifferentialEquations, LaTeXStrings, Symbolics

┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]

ArgumentError: Package DifferentialEquations not found in current path.
- Run import Pkg; Pkg.add("DifferentialEquations") to install the DifferentialEquations package.

Stacktrace:
 macro expansion
 macro expansion
@ ./lock.jl:223 [inlined]
 require(into::Module, mod::Symbol)
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


## 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

1. define the variables x_1 and x_2

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

3. build the function f_x with build_function

4. evaluate the first returned output from build_function

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

@variables x_1, x_2

LoadError: UndefVarError: @variables not defined
in expression starting at In:1

Stacktrace:
 top-level scope
@ :0
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


### 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]

UndefVarError: x_1 not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


### 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})

UndefVarError: build_function not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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

In this case, there are two functions:

1. 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]$$.

2. 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.0, 2.0)

UndefVarError: f_x not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

Plots.theme(:cooper)
x1_array = range(0, 1, length = 100)
x2_array = range(0, 1, length = 100)
Farray1 = [f_x(x1i, x2i) for x1i in x1_array for x2i in x2_array]
Farray2 = [f_x(x1i, x2i) 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, x1_array, x2_array, Farray2)

┌ Warning: :cooper is not a known theme, using :default
└ @ Plots /home/ryan/.julia/packages/Plots/GGa6i/src/themes.jl:10

UndefVarError: f_x not defined

Stacktrace:
 (::var"#1#4"{Float64})(x2i::Float64)
@ Main ./none:0
 iterate
@ ./generator.jl:47 [inlined]
 iterate
@ ./iterators.jl:1153 [inlined]
 iterate
@ ./iterators.jl:1147 [inlined]
 grow_to!(dest::Vector{Any}, itr::Base.Iterators.Flatten{Base.Generator{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, var"#2#3"}})
@ Base ./array.jl:860
 _collect
@ ./array.jl:764 [inlined]
 collect(itr::Base.Iterators.Flatten{Base.Generator{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, var"#2#3"}})
@ Base ./array.jl:712
 top-level scope
@ In:4
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


### 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])

UndefVarError: Symbolics not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


Next, build_function and eval the result.

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

UndefVarError: x_1 not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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

LoadError: UndefVarError: @variables not defined
in expression starting at In:1

Stacktrace:
 top-level scope
@ :0
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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

D = Differential(t)

UndefVarError: Differential not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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]

LoadError: UndefVarError: @variables not defined
in expression starting at In:1

Stacktrace:
 top-level scope
@ :0
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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

dz = D.(z)

UndefVarError: D not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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

g = [z*z, z+z+2*z, z^2-z^2]

UndefVarError: z not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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)

UndefVarError: D not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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))

UndefVarError: D not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


### 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(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,3,5,7], [2,4,6,8])

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

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


### Work-around for non-callable function variables¶

Here is my workaround:

1. create callable variables q and dq with the same length as z(t)

2. substitute the values into the Dg symbolic function and simplify the results

3. 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 to z

@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)),)))

LoadError: UndefVarError: @variables not defined
in expression starting at In:1

Stacktrace:
 top-level scope
@ :0
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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,3,5,7], [2,4,6,8])
# eval(Dg_func)([1,3,5,7], [2,4,6,8])

UndefVarError: build_function not defined

Stacktrace:
 top-level scope
@ In:1
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


## 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, t_space, sin.(t_space .+ phi), label = "z_$i") plot!(p, t_space, cos.(t_space .+ phi), label = "dz_$i/dt")
plot
end
p 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(z_space[i, :])
Dg_val[i, :] = Dg_func(z_space[i, :], dz_space[i, :])
end
for i in range(1, 3)
plot!(p, t_space, g_val[:, i], axis = nothing, label = "g_$i") plot!(p, t_space, Dg_val[:, i], axis = nothing, label = "dg_$i/dt")
end
p

UndefVarError: g_func not defined

Stacktrace:
 top-level scope
@ ./In:6
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)


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, t_space[1:ti], z_space[1:ti, i], axis = nothing, label = "", title="z", linecolor = color)
plot!(p, 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, t_space[1:ti], g_val[1:ti, i], axis = nothing, label = "", title = "g(z)",linecolor = color)
plot!(p, t_space[1:ti], Dg_val[1:ti, i], axis = nothing, label = "", title = "dg/dt", linecolor = color)
end
end

ArgumentError: Package ColorSchemes not found in current path.
- Run import Pkg; Pkg.add("ColorSchemes") to install the ColorSchemes package.

Stacktrace:
 macro expansion
 macro expansion
@ ./lock.jl:223 [inlined]
 require(into::Module, mod::Symbol)
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

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)

UndefVarError: g not defined

Stacktrace:
 top-level scope
@ :0
 eval
@ ./boot.jl:368 [inlined]
 include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

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$$.