2023-01-23 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, big(3.0) # integers, double precision, and single precision

(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

Arrays#

[1, 2, 3]

3-element Vector{Int64}:
 1
 2
 3

# explicit typing
Float64[1,2,3]

3-element Vector{Float64}:
 1.0
 2.0
 3.0

# promotion rules 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

# multi-dimensional array
A = [10 20 30; 40 50 60]

2×3 Matrix{Int64}:
 10  20  30
 40  50  60

A = np.array([[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:50000000

1:50000000

collect(1:5)

5-element Vector{Int64}:
 1
 2
 3
 4
 5

x = 0
for n in 1:50000000
    x += 1/n^2
end
@show x
x - pi^2/6 # trivia you can easily look up if needed

x = 1.6449340467988642

-2.0049362170482254e-8

# 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
(1 + 1.2e-16) - 1

2.220446049250313e-16

using Plots
default(linewidth=4)
plot(x -> (1 + x) - 1, xlims=(-1e-15, 1e-15),
xlabel="x", ylabel="y")
plot!(x -> x)
../_images/2023-01-23-functions_31_0.svg

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