#### Introduction

Julia is a functional language; that is, functions are first class objects that can be assigned to variables, passed as arguments, and returned from functions. This is a subtle, but powerful concept.

In this tutorial we will write a program that takes any second order polynomial and returns its roots. A simple problem for a human to solve, which becomes somewhat non-trivial when approached from a programming standpoint. We will approach the problem by creating a function that returns the derivative of another function and using that to calculate the coefficients of the polynomial, which can then be run through a quadratic formula.

#### Getting Started

Let's begin by creating a new file in Julia Studio. If this is your first time creating a file, you can do this through the file menu, by going to File > New File, or by clicking on the New File icon in the middle of the Welcome tab.

For this tutorial create a file called main.jl. This file will be where we begin execution of your program.

Now that we have a main file, let's create a new file called quadratic.jl, which will contain one of the methods that we will use to solve the problem. To create a file, use the File menu or the CTRL/CMD-N hotkey.

Now let's reference quadratic.jl from main.jl so it will be included in the execution of the program.

# include another file
# remember: only double quotes are valid in Julia


#### Background

We will approach this problem from the following angle:

Let's assume that we want to accept any second order polynomial as an argument, and we want to return its roots.

$f(x) = a{x}^2 + bx + c$

This can be easily represented with Julia:

f = (x) -> ax^2 + bx + c # a, b, and c are constants


We can then procede to find the coefficients of this polynomial by first, evaluating $f(0)$.

$c = f(0) = 0^2 * a + 0 * b + c$

It follows that, by taking the first derivative, we can obtain the value of b by evaluating $f’(0)$.

$f’ = \frac{df}{dx} = 2ax + b \\ b = f’(0)$

And the value of a can be obtained by evaluating $f(1)$ and subtracting a and b.

$f(1) = a + b + c \\ a = f(1) - b - c$

For now, let's assume that we have a method, which we will call derivative (more on this in a bit), that takes as an argument a function and returns another function, that represents the derivative of the function it is given. Read on to see how we implement this with Julia.

In quadratic.jl create a function called quadratic, which takes an argument y, which will be the quadratic equation we want to solve.

Unlike with MATLAB®, functions can be defined in any place, and called by any name. Although it is not required, for the sake of modularity, we usually follow the MATLAB convention of defining functions in a file of the same name.

# In quadratic.jl

# insert magic here
end


If you are new to functions in Julia you may have noticed that functions can be defined in multiple ways. Functions can be defined with the normal syntax (shown above), or as anonymous functions, for example:

# Examples of anonymous functions in Julia

square_stuff = (x) -> x^2

square_stuff = function(x)
return x^2
end

# The return statement is optional, but is generally a good idea for multi-line functions.
# By default the last line is returned if no return statement is found.


For more on functions, check out the documentation.

Let's write the code that actually gets the coefficients and plugs them into the good old quadratic formula, using a placeholder for our hypothetical derivative function for now.

function quadratic(f)
# Compute the first derivative of f
f1 = derivative(f)

# Get the y intercept (explictly passing a floating point number)
c = f(0.0)

# Get the y intercept of the first derivative
b = f1(0.0)

a = f(1.0) - b - c

# Our old friend the quadratic formula
# Notice how Julia lets you return multiple values
# Multiple values should be separated by a comma and
# are returned as a tuple
return (-b + sqrt(b^2 -4a*c)) / 2a, (-b - sqrt(b^2 -4a*c)) / 2a
end


That looks great, but what do we do if the roots are complex numbers? Not to worry, Julia has robust support for complex numbers. In this case, the easiest way to deal with complex roots is to simply add 0im to the determinant in our return values. That instructs Julia to treat the determinants as complex number types, and the sqrt function will handle them gracefully if they are negative. The Julia symbol for i, or $\sqrt{-1}$ is im, and can be prefixed with any numeric type to instantiate an imaginary number. In this example, we add 0im, which instantiates an imaginary number with magnitude 0. In our example, sqrt(b^2 -4a*c) becomes sqrt(b^2 -4a*c + 0im).

function quadratic(f)
# Compute the first derivative of f
f1 = derivative(f)

# Get the y intercept (explictly passing a floating point number)
c = f(0.0)

# Get the y intercept of the first derivative
b = f1(0.0)

a = f(1.0) - b - c

# Our old friend the quadratic formula
# Notice how Julia lets you return multiple values
# Multiple values should be separated by a comma and
# are returned as a tuple
return (-b + sqrt(b^2 -4a*c + 0im)) / 2a, (-b - sqrt(b^2 -4a*c + 0im)) / 2a
end


#### Numerical Differentiation

This is where Julia's functional nature helps it to shine. As you will see, the ability to receive functions as arguments and return them as results is very useful in this case.

You may recall the definition of a derivative:

$\frac{df}{dx} = \lim_{h\to0}\frac{f(x + h) - f(x)}{h}$

The question is, how do we express this in Julia, and how do we do it in such a way that we return a new function, rather than the derivative evaluated at a given point.

Let's start by creating a new file, which we will call derivative.jl. In this file, let's stub out a function.

# in derivative.jl

function derivative(f)
return f
end


This function satisfies the requirement of taking a function and returning the function, but at this point, it's not doing anything. Let's flesh this out a little more.

# in derivative.jl

function derivative(f)
# Look: you're returning a function of x

return function(x)  # yet another way to define functions in Julia
# pick a small value for h
# evaluate f at x + h
# evaluate f at x
# divide the difference by h
end
end

# derivative((x) -> x^2 + 2x + 3) returns a function that approximates (x) -> 2x + 2


Is it clear what we're doing there? When you call derivative, the function you pass in will be stored in a "closure". In other words, it will be stored so it can be called later. Derivative then returns another function, which takes x as an argument. The result of calling derivative will be a function, defined in terms of one variable, x. When that function gets called, it will use the stored function, f, to approximate the derivative at value, x. Let's fill in the gaps.

# in derivative.jl

function derivative(f)
return function(x)
# pick a small value for h
h = 0.00001

# evaluate f at x + h
f1 = f(x + h)

# evaluate f at x
f0 = f(x)

# divide the difference by h

return (f1 - f0) / h
end
end

# derivative((x) -> x^2 + 2x + 3) returns a function that approximates (x) -> 2x + 2


Theoretically, this is correct. However, notice the peculiar results based on your input after running it. This is due to limitations of floating point arithmetic, and a value of h which is ill-suited for the value of x being evaluated. See this Wikipedia article for a detailed discussion, but for our purposes, we will need to factor in machine epsilon to get precise answers. A good choice for a value of h to use is $\sqrt{\epsilon}x$, where is $\epsilon$ is machine epsilion. Conveniently, Julia provides a method, eps(Type) that gives us this value, right out of the box.

Also, due to limitations in floating point arithmetic, we'll have to do some gymnastics to avoid binary rounding errors.

# in derivative.jl

function derivative(f)
return function(x)
# pick a small value for h
h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x

# floating point arithmetic gymnastics
xph = x + h
dx = xph - x

# evaluate f at x + h
f1 = f(xph)

# evaluate f at x
f0 = f(x)

# divide the difference by h
return (f1 - f0) / dx
end
end

# derivative((x) -> x^2 + 2x + 3) returns a function that approximates (x) -> 2x + 2


#### Conclusion

Now we have all the code we need to make the program work. Let's wrap things up. In quadratic.jl we will need to load in the derivative method we wrote. The contents of the file should now be:

# Include our library
include("derivative.jl")