4) Introduction to programming in Julia#
Last time:
Introduction to PDEs
Today:
Using Julia
Errors
2.1. PlottingMachine Epsilon
1. Julia#
To recap from first class, Julia is a relatively new programming language. Think of it as MATLAB done right, open source, and fast. It’s nominally general-purpose, but mostly for numerical/scientific/statistical computing. There are great learning resources. We’ll introduce concepts and language features as we go.
It is supported in Jupyter Notebooks.
From your terminal, start a Julia session with
julia
To add the IJulia package to your environment, then type:
]add IJulia
go back to your Julia REPL using backspace and type
julia> using IJulia
and then
julia> notebook()
to run the notebook.
Let’s see it with live demos:
# The last line of a cell is output by default
x = 3
y = 4
4
println("$x + $y = $(x + y)") # string formatting/interpolation
3 + 4 = 7
4; # trailing semicolon suppresses output
Julia has several macros, invoked by the @ symbol. A useful one is the @show macro (which prints one or more expressions, and their results, to stdout, and returns the last result):
@show x + y # here the entire expression is shown
x * y # here only the return (computed) value is shown
x + y = 7
12
1.1. Numbers#
Tip
Check this very helpful Julia documentation page.
3, 3.0, 3.0f0, big(3.0) # integers, double precision, single precision, and convert to a maximum precision representation
(3, 3.0, 3.0f0, 3.0)
typeof(3), typeof(3.0), typeof(3.0f0), typeof(big(3.0))
(Int64, Float64, Float32, BigFloat)
# automatic promotion
@show 3 + 3.0
@show 3.0 + 3.0f0
@show 3 + 3.0f0;
3 + 3.0 = 6.0
3.0 + 3.0f0 = 6.0
3 + 3.0f0 = 6.0f0
# floating and integer division
@show 4 / 2
@show -3 ÷ 2; # type `\div` and press TAB
4 / 2 = 2.0
-3 ÷ 2 = -1
1.2 Arrays#
[1, 2, 3]
3-element Vector{Int64}:
1
2
3
# promotion rules for arrays are similar to arithmetic
[1,2,3.] + [1,2,3]
3-element Vector{Float64}:
2.0
4.0
6.0
x = [10., 20, 30]
x[2] # one-based indexing
20.0
x[2] = 3.5
3.5
Compare with the Python notation:
A = np.array([[10, 20, 30], [40, 50, 60]])
1.3 Functions#
# define a function called 'f' that takes three arguments, one of which has a default value
function f(x, y; z=3)
sqrt(x*x + y*y) + z
end
# define another function also called 'f' that only takes two arguments
function f(x, y)
sqrt(x*x + y*y)
end
f (generic function with 1 method)
What is the difference between a method and a function in Julia? Check this resource.
# if I invoke f this way, which one will the Julia compiler use?
f(3, 4, z=5)
10.0
# if instead I invoke f this way, which one will the Julia compiler use?
f(2, 3)
3.605551275463989
There are also compact “assignment form” function definitions:
g(x, y) = sqrt(x^2 + y^2)
g(3, 4)
5.0
In the assignment form above, the body of the function must be a single expression, although it can be a compound expression.
We have just seen an example of one of the most powerful features of Julia: dynamic (multiple) dispatch:
Functions can have multiple definitions as long each definition restricts the type of the parameters differently. It is the type of the parameters that define which “definition” (or “method” in Julia terminology) will be called.
This is the way Julia mimics polymorphism that some compiled languages (e.g., C++) have.
Anonymous functions are usually functions so short-lived that they do not need a name. This is done with an arrow notation. Example:
((x, y) -> sqrt(x^2 + y^2))(3, 4)
5.0
1.4 Ranges and loops#
# ranges
1:50000000
1:50000000
collect(1:5) # creates a vector from a range of values
5-element Vector{Int64}:
1
2
3
4
5
x = 0
# example of a for loop
for n in 1:50000000
x += 1/n^2
end
@show x
x - pi^2/6 # Basel problem (solved by Euler; generalized by Riemann's zeta function)
x = 1.6449340467988642
-2.0049362170482254e-8
# list comprehensions
sum([1/n^2 for n in 1:1000])
1.6439345666815612
2. Errors#
There are several sources of error in computation:
roundoff (a consequence of how computers represent numbers)
conditioning (a property of the problem we’re solving)
stability (a property of the algorithm we’re using)
2.1 Roundoff errors and Floating Point Arithmetic (FPA)#
If \(\circ\) is an operator that rounds a real number \(x \in \mathbb{R}\) to the nearest floating-point number \(\tilde{x}\):
Then the relative error introduced by this representation is
Rearranged:
What is a guaranteed upper bound on \(\delta x\)?
The IEEE (Institute of Electrical and Electronics Engineers) stadard guarantees that
where \(\varepsilon_M\) is machine precision.
Question: What is \(\varepsilon_M\) for single- and double-precision arithmetic?
eps(Float64)
2.220446049250313e-16
eps(Float32)
1.1920929f-7
It’s around \(10^{-7}\) for single-precision, and \(10^{-16}\) for double-precision floating point numbers.
What does this mean for the number of digits?
Why would anyone choose single precision?
x = 1e-16
@show 1+x # :-(
1 + x = 1.0
1.0
3. Plotting#
There are several plotting libraries in Julia. Plots.jl is a basic package for 1D plots and simple planar maps. Makie.jl is a much more sophisticated package with richer 2D/3D visualizations, dynamic sliders, animations, etc.
We will use Plots.jl in class.
Before being able to use it, you have to add it to your environment.
using Pkg
Pkg.add("Plots") # adds the Plots.jl package to your environment
using Plots # once added, you can use it
default(linewidth=4, legendfontsize=12, xtickfontsize=12, ytickfontsize=12) # sets some defaults for embelishments
Resolving package versions...
No Changes to `~/.julia/environments/v1.10/Project.toml`
No Changes to `~/.julia/environments/v1.10/Manifest.toml`
┌ Warning: Circular dependency detected.
│ Precompilation will be skipped for dependencies in this cycle:
│ ┌ Symbolics → SymbolicsForwardDiffExt
│ └─ Symbolics → SymbolicsPreallocationToolsExt
└ @ Pkg.API.Precompilation ~/julia-1.10.10/share/julia/stdlib/v1.10/Pkg/src/precompilation.jl:583
x = range(0, 10, length=100)
y = sin.(x)
plot(x, y, label = "f(x)=sin(x)", xlabel = "x", ylabel = "y")
