2022-01-12 Functions

  • Using Julia

  • Plotting

  • Intro to floating point

  • Summing series

  • Relative vs absolute errors

Julia

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.

# The last line of a cell is output by default
x = 3
y = 4
4
println("$x + $y = $(x + y)") # string formatting/interpolation
4;  # trailing semicolon suppresses output
3 + 4 = 7
@show x + y
x * y
x + y = 7
12

Numbers

3, 3.0, 3.0f0 # integers, double precision, and single precision
(3, 3.0, 3.0f0)
typeof(3), typeof(3.0), typeof(3.0f0)
(Int64, Float64, Float32)
# 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 3 / 2
@show -3 ÷ 2; # type `\div` and press TAB
3 / 2 = 1.5
-3 ÷ 2 = -1

Arrays

[1,2,3]
3-element Vector{Int64}:
 1
 2
 3
# explcit typing
Float64[1,2,3]
3-element Vector{Float64}:
 1.0
 2.0
 3.0
# promotion rules similar to arithmetic
[1,2,3.]
3-element Vector{Float64}:
 1.0
 2.0
 3.0
x = [10, 20, 30]
x[2] # one-based indexing
20
#x[2] = 3.5
# multi-dimensional array
A = [10 20 30; 40 50 60]
2×3 Matrix{Int64}:
 10  20  30
 40  50  60

Functions

function f(x, y; z=3)
    sqrt(x*x + y*y) + z
end
f (generic function with 1 method)
f(3, 4, z=5)
10.0
g(x, y) = sqrt(x^2 + y^2)
g(3, 4)
5.0
((x, y) -> sqrt(x^2 + y^2))(3, 4)
5.0

Loops

# ranges
1:5
1:5
collect(1:5)
5-element Vector{Int64}:
 1
 2
 3
 4
 5
x = 0
for n in 1:5000
    x += 1/n^2
end
@show x
x - pi^2/6 # trivia you can easily look up if needed
x = 1.6447340868469014
-0.00019998000132503968
# list comprehensions
sum([1/n^2 for n in 1:1000])
1.6439345666815612

Poll: What is floating point arithmetic?

  1. fuzzy arithmetic

  2. exact arithmetic, correctly rounded

  3. the primary focus of numerical analysis

0.1 + 0.2
0.30000000000000004
using Plots
default(linewidth=4)
plot(x -> (1 + x) - 1, xlims=(-2e-16, 2e-16))
plot!(x -> x)
../_images/2022-01-12-functions_30_0.svg
┌ Warning: No strict ticks found
└ @ PlotUtils /home/jed/.julia/packages/PlotUtils/xekml/src/ticks.jl:333
┌ Warning: No strict ticks found
└ @ PlotUtils /home/jed/.julia/packages/PlotUtils/xekml/src/ticks.jl:333
┌ Warning: No strict ticks found
└ @ PlotUtils /home/jed/.julia/packages/PlotUtils/xekml/src/ticks.jl:333
┌ Warning: No strict ticks found
└ @ PlotUtils /home/jed/.julia/packages/PlotUtils/xekml/src/ticks.jl:333

Machine epsilon

We approximate real numbers with floating point arithmetic, which can only represent discrete values. In particular, there exists a largest number, which we call \(\epsilon_{\text{machine}}\), such that

\[ 1 \oplus x = 1 \quad \text{for all}\ \lvert x \rvert < \epsilon_{\text{machine}}.\]

The notation \(\oplus, \ominus, \odot, \oslash\) represent the elementary operation carried out in floating point arithmetic.

eps = 1
while 1 + eps != 1
    eps = eps / 2
end
eps
1.1102230246251565e-16
eps = 1.f0
while 1 + eps != 1
    eps = eps / 2
end
eps
5.9604645f-8

Beating exp

\[e^x = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \dotsb\]

Suppose we want to compute \(f(x) = e^x - 1\) for small values of \(x\).

f1(x) = exp(x) - 1
y1 = f1(1e-8)
9.99999993922529e-9
f2(x) = x + x^2/2 + x^3/6
y2 = f2(1e-8)
1.000000005e-8

Which answer is more accurate?

@show (y1 - y2)        # Absolute difference
@show (y1 - y2) / y2;  # Relative difference
y1 - y2 = -1.1077470910720506e-16
(y1 - y2) / y2 = -1.1077470855333152e-8
0 / 0
NaN