5) Functions#
Last time:
Introduction to Fortran
Today:
Using Julia
Introduction to floating point arithmetic
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.
# 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
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
# 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
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
Argument Passing Behavior#
Julia function arguments follow a convention sometimes called “pass-by-sharing”, which means that values are not copied when they are passed to functions. Function arguments themselves act as new variable bindings (new “names” that can refer to values), much like assignments argument_name = argument_value
, so that the objects they refer to are identical to the passed values. Modifications to mutable values (such as Array
s) made within a function will be visible to the caller. (This is the same behavior found in Scheme, most Lisps, Python, Ruby and Perl, among other dynamic languages.)
As a common convention in Julia (not a syntactic requirement), a function that modifies any of its arguments would typically be named with a !
, e.g., f!(x, y)
rather than f(x, y)
, as a visual reminder at the call site that at least one of the arguments (often the first one) is being mutated.
To read more about Julia functions, check the Julia documentation.
Argument-type declarations#
You can declare the types of function arguments by appending ::TypeName
to the argument name. Example:
function f(x::Float64, y::Float64)
sqrt(x*x + y*y)
end
f (generic function with 2 methods)
Note: Argument-type declarations normally have no impact on performance: regardless of what argument types (if any) are declared, Julia compiles a specialized version of the function for the actual argument types passed by the caller.
The most common reasons to declare argument types in Julia are, instead: dispatch, correctness, and clarity.
1.4 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 # 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. Floating point arithmetic#
0.1 + 0.2
(1 + 1.2e-16) - 1
2.220446049250313e-16
using Pkg
Pkg.add("Plots")
using Plots
default(linewidth=4)
plot(x -> (1 + x) - 1, xlims=(-1e-15, 1e-15),
xlabel="x", ylabel="y")
plot!(x -> x)
Resolving package versions...
No Changes to `~/Dropbox/SDSU/Teaching/Comp526/fall25/Project.toml`
No Changes to `~/Dropbox/SDSU/Teaching/Comp526/fall25/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.9/share/julia/stdlib/v1.10/Pkg/src/precompilation.jl:583
3. 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
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
3.1 Approximating exp
#
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