The Julia programming language

Numerical Analysis

Urbain Vaes

NYU Paris

Course organization (1/3)

Organisation of each class

Each class is divided into two parts:

  • a theoretical lecture, usually with accompanying slides;

  • a recitation comprising pen-and-paper exercises and implementation tasks.

\(~\)

Evaluation methods

  • Jupyter notebooks to hand in, roughly one per chapter (25%)

  • Weekly 15-minute quizzes, based on multiple-choice questions (25%)

  • Midterm exam with theoretical and practical questions, without AI (Bonus for exam)

  • Final exam with theoretical and practical questions, without AI (50%): \[ E = \max \{ E_{final}, 0.4 \times E_{midterm} + 0.6 \times E_{final} \} \]

Course organization (2/3)

Detailed program

  • Chapter 0 : Introduction to Julia
  • Chapter 1 : Floating point arithmetic
  • Chapter 2 : Interpolation and integration
  • Chapter 3 : Linear systems
  • Chapter 4 : Eigenvalue problems
  • Chapter 5 : Nonlinear systems
  • Chapter 6 : Ordinary differential equations


Important remarks

  • The first notebook is essential for the rest of the course.

  • The course advances quickly and requires regular work.

  • Unjustified absences will be handled according to NYU policy.

Course organization (3/3)

Sources of error in computational science

  • Modeling error: discrepancy between the mathematical model and reality
  • Data error: uncertainty in the data of the problem
  • Discretization error: error when turning mathematical equations into solvable finite-dimensional problems
  • Discrete solver error: error introduced by the numerical method used to solve the discretized equations
  • Round-off errors: limited accuracy of computer arithmetics

The latter three errors, particularly discrete solver errors, will be the focus of this course.

⚠️ Timetable adjustments

  • The class of Wed 3 Sep is postponed to the morning of Friday 12 Sep

  • The class of Mon 15 Sep is postponed to the morning of Friday 19 Sep

Julia, a good language for scientific computing?

Python Matlab Julia C
Libre? ✔️ ❌ ✔️ ✔️
Fast? ❔ ❔ ✔️ ✔️
Easy? ✔️ ✔️ ✔️ ❌
Math-friendly? ❔ ✔️ ✔️ ❌
Ecosystem? ✔️ ❔ ❔ ✔️

  • Julia is easy to learn, read and maintain

  • Julia is fast: “As simple as Python with the speed of C++”

  • Julia is well-suited for numerical mathematics

  • Julia is libre software

A bit of history

  • Created in 2009 by J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman

  • Excerpt from the Julia Manifesto:

We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled.

Adoption of Julia is rapidly increasing

Performance benchmark

Julia is a great fit for numerical mathematics

  • Elegant syntax to define functions
f(x) = cos(x) + sin(x)
f(1)
1.3817732906760363
  • The syntax for anonymous functions is equally elegant:
import Plots
Plots.plot(x -> cos(x) + sin(x))
  • Use of unicode characters:
ω = 1  # ⚠️ Even emojis work!
f(t) = sin(ω*t)
∇f(t) = ω * cos(ω*t)
f²(t) = f(t) * f(t)
  • Compiled just in time (for loops allowed!):
@elapsed sum(1/n^2 for n ∈ 1:10^8)
0.131086146

In comparison, Python takes about 10 seconds…

  • Contains most math functions out of the box:
π, Float64(ℯ), randn(), √2
(π, 2.718281828459045, -0.924364916359138, 1.4142135623730951)
  • Concise syntax, MATLAB-like syntax for matrices:
M = [1 2 3; 4 5 6]
2×3 Matrix{Int64}:
 1  2  3
 4  5  6
M[:, [1, 3]]
2×2 Matrix{Int64}:
 1  3
 4  6
  • Great syntax for array broadcasting:
broadcast(gcd, 6, [9, 8, 7])
gcd.(6, [9 8 7])
1×3 Matrix{Int64}:
 3  2  1
1 .+ [1; 2; 3] .+ [1 1 1; 1 1 1; 1 1 1]
3×3 Matrix{Int64}:
 3  3  3
 4  4  4
 5  5  5

  • Dynamic types by default…
x = 1 ; y = π ; z = 1//2 ; t = x + y + z
4.641592653589793
for v ∈ (x, y, z, t)
    println("$v is a $(typeof(v))")
end
1 is a Int64
π is a Irrational{:π}
1//2 is a Rational{Int64}
4.641592653589793 is a Float64
  • … but static typing allowed
f(x) = 2x
f(x::Float64) = 3x
@show f(1)
@show f(1.) ;
f(1) = 2
f(1.0) = 3.0
  • Array comprehension (for within array)
A = [63/(i+2j) for i ∈ 1:3, j ∈ 1:3]
3×3 Matrix{Float64}:
 21.0   12.6  9.0
 15.75  10.5  7.875
 12.6    9.0  7.0
  • Built-in linear solver
x = [1.2, 3.4, 5.6]
b = A*x
A\b
3-element Vector{Float64}:
 1.2000000000000222
 3.3999999999999084
 5.600000000000077
@show A\b == x
@show A\b ≈ x ;
A \ b == x = false
A \ b ≈ x = true
  • Array indices start at 1
a = [1, 2, 3]
@show a[1]
@show a[2]
@show a[3]
a[1] = 1
a[2] = 2
a[3] = 3
3
  • New binary operators
⊕(α,β) = α + β
5 ⊕ 7
12

Julia’s programming paradigm: multiple dispatch

Wikipedia definition: paradigm in which a function or method can be dynamically dispatched based on the run-time type.

  • Example:
f(x::Int) = "Int $x"
f(x::Float64) = "Float: $x"
f(x::Any) = "Generic fallback"
f(1.2)
"Float: 1.2"
  • Makes it simple to extend existing packages:
import Base.-
function -(s1::String, s2::String)
    return s1[1:length(s1)-length(s2)]
end
"julia" - "ia"
"jul"
  • Python, in comparison, uses single dispach.

The REPL and Julia programs

  • REPL = read–eval–print loop = “command line”

  • To access the REPL,

    • either type julia in a terminal window;

    • or launch it from VSCode by typing Shift+Ctrl+P then Julia: Start REPL

  • Useful shortcuts in the REPL:

    • ] to access package mode

    • ; to access shell mode

    • ? to access help mode

    • ⟵ to return to normal mode

    • Tab for automatic completion

    • up and down arrows to navigate history

  • To write a julia program

    • either write a .jl file, which can be run using Shift+Enter

    • Or write a Jupyter notebook with .ipynb extension.

Package management

To install a new package

  • Enter package mode with ;

  • Type add followed by the package name, e.g.

    add Plots

  • You only need to do this once

To use an installed package

  • Either use import

    import Plots
Plots.plot(cos)

  • or use using (no need to prepend package name later)

    using Plots
plot(cos)

Additional references