Numerical methods in Julia tutorial

Introduction

Welcome to JuliaTutorial.jl! The aim of this page is to get you quickly started with the Julia programming language.

Julia is perhaps the most productive gateway to the world of numerical methods, and I hope this can be a self-contained resource to get you started with Julia, its ecosystem, and exploring numerical methods with a toy example of making a linear solver using LU decomposition.

The LU decomposition tutorial was largely adapted from a section in Kyle Novak's book Numerical Methods for Scientific Computing (Chapter 2: Direct Methods for Linear Systems), so I highly recommend that as a supplementary resource.

Installation

Julia can be installed from its official download page.

It can also be installed via homebrew (on macs) by pasting the following command into a terminal. This command will install the current stable release (currently 1.8.1).

brew install --cask julia

Simply type julia in a terminal to open a Julia REPL (read-evaluate-print-loop) session. You can use the REPL to evaluate simple Julia expression, manage Julia packages, run shell commands, and access documentation. Additionally, you can run a Julia file in you active directory by typing julia my_julia_program.jl.

Let's open a julia REPL and type print("hello world").

~ % julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.8.1 (2022-09-06)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> print("hello world")
hello world

The REPL is a great interface for using Julia quickly (to make a fast plot, overview your packages, or do a small calculation), but we will need a robust IDE (Integrated Development Environment) to write more complex code. There are many IDE's for various programming languages including Julia, but the easiest and most fully featured one for Julia is VS Code which we will install in the next section.

Using Julia in VS Code

VS Code can be installed from its official download page.

It can also be installed via homebrew (on macs).

brew install --cask visual-studio-code

You can then open VS Code by typing code into a terminal.

1. The Julia extension can be install by clicking the following extension button in the upper left of VS Code.

1. Search for "julia", and then click the blue install button.

Next, let's make a .jl file to edit. Click new file.

And select Julia File to initialize a blank Julia file.

Finally, let's open up a REPL inside of VS Code (called the integrated terminal) using the default hotkey ctrl + J + O.

Congratulations! We are now set-up for working with Julia.

Writing simple Julia programs

Type 1+1 into the editor window in VS Code and press shift + enter to execute the line your cursor is on (or execute multiple lines if you have them highlighted). By default, it will print an "inline result" next to the code as well as in the REPL session below your editor.

This lends itself to an interactive workflow much like a jupyter notebook. This is the most basic workflow and suitable for writing scripts and exploring data, but if you are doing package development with unit tests and more complicated environments you will want to look at Revise.jl.

The Julia package manager

Julia packages such as DifferentialEquations.jl or Plots.jl are collections of functions which can be imported to provide additional functionality to your Julia program. They are equivalent to libraries in other languages, for instance, python's numpy and Julia's LinearAlgebra.jl provide similar functionality.

The Julia package manager, Pkg.jl, is itself a package, and it lets us import packages from various locations on the internet and manage them. It is essentially python's pyp and conda/venv all in one.

We can access Pkg.jl in two ways.

1. This is the less common way, but it can be accessed like any Julia package by typing using Pkg in either your editor window (with shift+enter) or the REPL. You can then import other packages (here Plots.jl) by typing Pkg.add("Plots.jl"). This will download Plots.jl from its official github repo to your ~/.julia/packages folder and updates your project.toml file in .julia/environments/1.8 which contains the default v1.8 environment. You can then access Plots.jl functions at any time by running using Plots.jl in Julia.

2. The more common way is using the built in REPL interface. Type ] into your REPL window. The julia> prompt will change into (@v1.8) pkg> indicating it now accepts Pkg.jl functions (use backspace to return to julia>). We can type add Plots to do the same thing we did above. In addition to add you can also remove/rm packages, list all your current packages with status/st, or run unit tests for the package test Plots.

Plotting

Creating plots is perhaps the most fundamental thing you can do in scientific computing and it is also probably the quickest and most fun way to learn Julia.

To make a plot in Julia, run the following.

using Plots
plot(sin)

We can also call plot on an array.

plot([1,2,4,2,5])

Or put one plot on top of another by calling plot!. In Julia, it is convention to put a "!" after a function which modifies data. In this case, plot!([1,2,4,2,5]) writes over our previous plot(sin).

plot(sin)
plot!([1,2,4,2,5])

Understanding the plot function call

Let's take a minute to unpack the function call plot(sin). The function plot is explicitly exported from the "Plots.jl" package, so we do not need to write Plots.plot to specify it (though we would if we if we were working with several different plotting packages at once).

This function consists of many "sub-functions" known as "methods" which are specialized to different arguments for the function. When you call the function plot(sin), Julia sees that you called plot on a single argument sin of type Function (functions are treated as just another data type in Julia such as Int64, FLoat64, etc). It can then use a method of plot which calls the function sin on a few values (here -5 to 5) which serve as the values on the x-axis.

In contrast, plot([1,2,4,2,5]) uses a slightly different method where it generates (x,y) coordinates based off of the index of the value in the array and the actual value of the array respectively.

It is worth pointing out here that unlike python and many other programming languages, Julia has 1-based arrays meaning the first index is 1 not 0 which is why you see the array plotted as 1 to 5 on the x-axis.

Making animations with a few important language features

It is also possible to make simple animations in Plots.jl.

Let's make a sine wave which translates itself to the right.

@gif for t in 0:0.5:6
    plot(x->sin(x-t))
end

There is a lot of syntax to unpack here, but I think it introduces many important language features of Julia.

1. @gif is a "macro"
2. for t in 0:0.5:6...end is a "for loop"
3. x->sin(x-t) is an "anonymous function" or "lambda"

Macros

Macros in Julia can be thought of as generalizations of functions.

A function takes arguments which are various data types for instance sin takes in types of Number such as Int64, Float64, etc and gives you a value you would associate with sin from its mathematical definition.

But instead of taking normal data types like Int64, a macro takes a piece of Julia code as an argument. This is an example of what is known as "metaprogramming". In this case, @gif uses our for loop with a plot call inside as a recipe to make an animation.

You will likely rarely write these, but you will encounter and use them everywhere in the Julia ecosystem.

The most useful macro in Base Julia is @time which will time the code after it and write it to stdout (usually just the REPL). You can read the documentation for a macro by using ?@time in a REPL.

@time sin(pi/2)

  0.000001 seconds

julia> @btime sin(pi/2)
1.167 ns (0 allocations: 0 bytes)

julia> @time @benchmark plot(sin)
10.841423 seconds (247.36 M allocations: 5.504 GiB, 10.66% gc time, 0.31% compilation time)
BenchmarkTools.Trial: 5812 samples with 1 evaluation.
Range (min … max):  786.833 μs …   6.679 ms  ┊ GC (min … max): 0.00% … 80.84%
Time  (median):     805.958 μs               ┊ GC (median):    0.00%
Time  (mean ± σ):   858.419 μs ± 416.177 μs  ┊ GC (mean ± σ):  3.84% ±  6.87%

██▆▄▃▃▂▂▁▁ ▁▁                                                 ▁
██████████████▇▇▇▇▆▆▆▅▅▅▄▄▅▃▅▁▃▄▄▃▁▁▃▁▃▁▃▁▄▁▁▃▃▁▁▁▃▁▁▁▁▁▁▁▃▁▃ █
787 μs        Histogram: log(frequency) by time        1.6 ms <

Memory estimate: 497.41 KiB, allocs estimate: 21325.

For loops

Most programming languages have constructions called "for loops" which allow execution of code repeatedly.

In Julia the typical syntax is for i in Iterators...end. Note "white space" is not as important in Julia as it is in python for i in 1:5 end is a valid loop in Julia. This is the fastest and most flexible way of making a for loop and is therefore the most commonly used for writing numerical methods.

i is the index variable (can be any other letter(s)) and takes on the value of every element in "Iterator". Iterators are any data structure with many elements and an a ordering such as arrays, strings, and ranges. Here 1:5 is of type Range. We could also specify 1:2:5 to count by 2's or range(0,2pi,length=5) to go from 0 to 2pi in 5 steps.

for i in 1:5
    print(i)
end

12345

julia> [sin(i) for i in 1:3]
3-element Vector{Float64}:
0.8414709848078965
0.9092974268256817
0.1411200080598672

julia> (i for i in 1:3)
Base.Generator{UnitRange{Int64}, typeof(identity)}(identity, 1:3)
julia>sum(i for i in 1:3)
6

    sin.(1:3)
3-element Vector{Float64}:
0.8414709848078965
0.9092974268256817
0.1411200080598672

julia> map(sin, 1:3)
3-element Vector{Float64}:
0.8414709848078965
0.9092974268256817
0.1411200080598672

Anonymous functions

In many cases, we will want to define small functions that are essentially only used once in a program. Rather than clutter everything up with their definitions, we can use anonymous functions also sometimes called "lambdas".

These are written like x->sin(x) or (x,t)->sin(x-t) for multiple arguments. Compare with python lambda x: sin(x).

By putting them in parentheses, we can call them like a normal function.

julia> (x->sin(x))(pi/2) == sin(pi/2)
true

f(x,y) = plot(sin(cos(atan(y/x)+2)^2)

n = cos(atan(y/x)+2)^2
b = sin(n)
p = plot(b)

f(x,y) = atan(y/x) |> _->cos(_+2)^2 |> sin |> plot

The Julia ecosystem

The Julia ecosystem has ~8000 packages in the official registry. These can be viewed and searched on JuliaHub, and any one can be imported through Pkj.jl add.

I think one good informal place to start to get a feel for the scope of the current Julia community and its ecosystem is by looking at the playlist of presentations for JuliaCon 2022. We can see a strong focus primarily in applied mathematics, GPU and distributed computing, physics, economics, teaching, and statistics/data science.

Here is also a lineup of some notable packages that physical chemistry people might be interested in.

Scientific Computing

• SciML: This is a huge, unified ecosystem for doing scientific machine learning. Its documentation splash page is also a good place to discover scientific computing packages both in and outside its framework.

• DifferentialEquations.jl: The most notable member of the SciML ecosystem and probably the most comprehensive framework for numerically solving ordinary differential equations out there.

• Flux.jl: This is the main framework for machine learning in Julia and provides similar functionality to python's PyTorch. Compared to python machine learning frameworks, it is extremely lightweight and explicit. In fact, Flux.jl's documentation describes how you could have built Flux.jl yourself by simply making convenience functions for regressions using automatic differentiation.

• Zygote.jl: Julia's main automatic differentiation package and is what powers Flux.jl. It has a similar lightweight/explict design philosophy. The book Numerical Methods for Scientific Computing has a great tutorial on how to build an automatic differentiation tool in Julia similar to Zygote.

• Pluto.jl: Fills a similar role to Jupyter notebooks (actually you can use Julia in Jupyter as well).Has many advantages over Jupyter notebooks notably being just a normal .jl file as opposed to having both .py and .ipynb file types and being reactive - it updates all affected cells after changing a variable. Also, Grant Sanderson (3Blue1Brown) is a fan if you need any more reason to check it out.

Plotting

• Plots.jl: The jack-of-all-trades standard plotting package. Has many backends and a simple syntax to make quick and dirty plots.
• Makie.jl: An extremely flexible GPU powered plotting package. Can make nice statistical plots as well as interactive interfaces for "plots" which are borderline simple video games.
• Gadfly.jl: Specialized for making really clean, composable statistical plots.

Computer Algebra Systems

• MathLink.jl: Wolfram Mathematica has a completely free and lightly licensed version known as WolframEngine (can easily install it with brew). You can use it independently, but it is a really excellent interface in Julia that can bridge the gap between computational and analytic methods though at some performance cost. Combine with Pluto.jl for a pretty nice and very free version of a Mathematica notebook.
• Symbolics.jl: The pure Julia equivalent of Mathematica. It is not as user friendly or as flexible as Mathematica particularly with symbolic integration, but it is much faster and supports auto-differentiation.
• MetaTheory.jl: The backend for Symbolics.jl. You can define generic relationships between data structure in an abstract algebra-y way such as introducing an identity operation or a distributive property. It can then use Julia's metaprogramming interface and a concept known as E-graph saturation to simplify expressions of these structures.

Quantum

• QuantumOptics.jl: This is the biggest quantum package in the Julia ecosystem and similar to python's QuTip.

• OpenQuantumTools.jl: Contains performant solvers, mostly using DifferentialEquations.jl for working with open quantum systems.

a((u,p),(v,q)) =
    ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )*dΩ

l((v,q)) =
    ∫( v⋅f + q*g )*dΩ +
    ∫( v⋅(n_Γ⋅∇u) - (n_Γ⋅v)*p )*dΓ

    julia> π ≈ pi
    true

    julia> for 😄 ∈ [:α², :🐈,:🐕,:🐘, :ℶₐ] ∩ [:🐈,:🐕,:🐘] print(😄) end
    🐈🐕🐘