Lecture 38: Introduction to Julia¶
Julia is a programming language which aims to combine the high productivity of scripting languages with the high performance of compiled languages. The paper [BEKS17] describes the language.
Fig. 93 The logo of Julia, from the Software Engineering Daily web site.¶
Getting Started¶
In the cloud, at CoCalc: https://cocalc.com you can open a Jupyter notebook with a Julia kernel.
To install Julia on your own computer, follow the instructions
at <https://julialang.org>. To execute Julia programs in a Jupyter notebook,
do using IJulia followed by notebook() in the REPL of Julia,
that is: when running in a Terminal or PowerShell window.
With using one can install Plots, SymPy and Symbolics.
The SymPy.jl is a wrapper to the Python package SymPy,
whereas Symbolics [GMCGSER22] is a computer algebra system
written in Julia.
To demonstrate symbolic computation with Julia, consider the problem of computing the slope of the tangent line at a point on a circle.
Implicit Differentiation with SymPy.jl¶
In the equation of the circle \(x^2 + y^2 - 1 = 0\), the variable \(x\) is the independent variable and the variable \(y\) depends on \(x\). We will obtain the slope via implicit differentiation, as \(\displaystyle \frac{d}{d~\!x} y(x)\).
using SymPy
If SymPy is not yet installed, then you will be prompted to install.
The declarations of the variables result in dy to be assigned the
unknown slope \(\displaystyle \frac{d}{d~\!x} y(x)\).
x = Sym("x")
y = SymFunction("y")
dy = diff(y(x), x)
Observe that y is declared as a symbolic function,
in any symbol. If we want to use y when we define
the equation, then we must use y(x) to indicate that
the function y depends on x.
c = x^2 + y(x)^2 - 1
dc = diff(c, x)
dydx = solve(dc, dy)
The solve returns an Array which starts in Julia with index one.
The slope is then obtained as
symbolicslope = dydx[1]
and equals -x/y(x).
Let us now take a random point on the circle, generated by a random angle, using the polar representation.
angle = 2*pi*rand()
(a, b) = (cos(angle), sin(angle))
sd = Dict(x => a, y(x) => b)
The coordinates of the point (a, b) are stored in a dictionary,
which is a convenient data structure for substitution.
Observe that, as Julia is geared towards numerical computing,
the pi is a 64-bit floating-point approximation for \(\pi\).
The numerical value for the slope is then obtained via substitution
as follows:
slope = symbolicslope.subs(sd)
The numerical value in slope serves as a reference to check
the alternative computation of the slope via the Taylor series.
The series of SymPy works for one variable only,
but it works well on purely symbolic expressions, so we can compute
the expression for the tangent line by running series twice.
Because we continue the session, we declare new symbols X and Y
for the variables in the circle equation and A and B for the
coordinates of the point.
X, Y, A, B = Sym("X, Y, A, B")
nc = X^2 + Y^2 - 1
cx = SymPy.series(nc, X, A, 2)
In cx is the Taylor series expansion of the expression
stored in nc about X = A, of order 2.
To convert the cx into a symbolic expression on which
we can apply series again, we select the coefficient
to make the expression px (a polynomial as truncated series).
px = cx.coeff(X, 0) + cx.coeff(X - A, 1)*(X - A)
cy = series(px, Y, B, 2)
cY = cy.subs(Dict(A => a, B => b)
Then the slope is obtained as
- cY.coeff(X, 1)/cY.coeff(Y, 1)
which agrees with the value earlier assigned to slope.
Implicit Differentiation with Symbolics.jl¶
We will now compute the value of the slope a third time.
using Symbolics.jl
@variables x, y(x)
Observe that with Symbolics (opposite to what we did with SymPy),
we declare y(x) so y is explicitly defined as depending on y.
This implies that typing y suffices and, unlike with SymPy,
we do not need to type y(x) when we define the equation of the circle.
In order not to confuse with the Differential of SymPy,
which we used earlier in this notebook,
we declare that we use the Differential of Symbolics.
D = Symbolics.Differential(x)
dy = D(y)
As before, dy
is \(\displaystyle \frac{d}{d~\!x} y(x)\).
c = x^2 + y^2 - 1
dc = D(c)
The differential operator D does not expand automatically,
we have to ask for it.
eq = expand_derivatives(dc)
SBslope = Symbolics.solve_for(eq,dx)
In SBslope is the expression -x/y(x)
so the symbolic part of this problem is thus indeed solved.
What is left is the substitution, again using a dictionary,
mapping the coordinates (a, b) of the random point
to the symbols x and y.
nsd = Dict(x => a, y => b)
substitute(SBslope, nsd)
and we see again the value that was earlier assigned
to slope.
References
J. Bezanson, A. Edelman, S. Karpinski, V.B. Shah. A fresh approach to numerical computing. SIAM Review 59(1): 65-98, 2017.
S. Gowda, Y. Ma, A. Cheli, M. Gwozdz, V.B. Shah, A. Edelman, and C. Rackauckas. High-performance symbolic-numerics via multiple dispatch, arxiv:2105.03949v3 [cs.CL] 5 Feb 2022.