# Day_07 - Julia’s Symbolics: Computer Algebra System

## Contents

# 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]
└ @ Base loading.jl:1662
```

```
ArgumentError: Package DifferentialEquations not found in current path.
- Run `import Pkg; Pkg.add("DifferentialEquations")` to install the DifferentialEquations package.
Stacktrace:
[1] macro expansion
@ ./loading.jl:1163 [inlined]
[2] macro expansion
@ ./lock.jl:223 [inlined]
[3] require(into::Module, mod::Symbol)
@ Base ./loading.jl:1144
[4] eval
@ ./boot.jl:368 [inlined]
[5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
LoadError: UndefVarError: @variables not defined
in expression starting at In[2]:1
Stacktrace:
[1] top-level scope
@ :0
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

### 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:
[1] top-level scope
@ In[3]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

### 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 withwhereandwhenthe function is defined. This method using`eval`

creates a newworldwithin 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:
[1] top-level scope
@ In[4]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
UndefVarError: f_x not defined
Stacktrace:
[1] top-level scope
@ In[5]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
┌ 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:
[1] (::var"#1#4"{Float64})(x2i::Float64)
@ Main ./none:0
[2] iterate
@ ./generator.jl:47 [inlined]
[3] iterate
@ ./iterators.jl:1153 [inlined]
[4] iterate
@ ./iterators.jl:1147 [inlined]
[5] 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
[6] _collect
@ ./array.jl:764 [inlined]
[7] collect(itr::Base.Iterators.Flatten{Base.Generator{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, var"#2#3"}})
@ Base ./array.jl:712
[8] top-level scope
@ In[6]:4
[9] eval
@ ./boot.jl:368 [inlined]
[10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

### 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:
[1] top-level scope
@ In[7]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

Next, `build_function`

and `eval`

the result.

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

```
UndefVarError: x_1 not defined
Stacktrace:
[1] top-level scope
@ In[8]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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[9]:1
Stacktrace:
[1] top-level scope
@ :0
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
D = Differential(t)
```

```
UndefVarError: Differential not defined
Stacktrace:
[1] top-level scope
@ In[10]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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[11]:1
Stacktrace:
[1] top-level scope
@ :0
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
dz = D.(z)
```

```
UndefVarError: D not defined
Stacktrace:
[1] top-level scope
@ In[12]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
UndefVarError: z not defined
Stacktrace:
[1] top-level scope
@ In[13]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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:
[1] top-level scope
@ In[14]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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 afeaturenot abug.

```
Dg = expand_derivatives.(D.(g))
```

```
UndefVarError: D not defined
Stacktrace:
[1] top-level scope
@ In[15]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

### 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 resultsbuild 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)),)))
```

```
LoadError: UndefVarError: @variables not defined
in expression starting at In[16]:1
Stacktrace:
[1] top-level scope
@ :0
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
UndefVarError: build_function not defined
Stacktrace:
[1] top-level scope
@ In[17]:1
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

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

```
UndefVarError: g_func not defined
Stacktrace:
[1] top-level scope
@ ./In[20]:6
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

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

```
ArgumentError: Package ColorSchemes not found in current path.
- Run `import Pkg; Pkg.add("ColorSchemes")` to install the ColorSchemes package.
Stacktrace:
[1] macro expansion
@ ./loading.jl:1163 [inlined]
[2] macro expansion
@ ./lock.jl:223 [inlined]
[3] require(into::Module, mod::Symbol)
@ Base ./loading.jl:1144
[4] eval
@ ./boot.jl:368 [inlined]
[5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

```
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:
[1] top-level scope
@ :0
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base ./loading.jl:1428
```

**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\).