The julia language is a new language and as such, certain design decisions are still being made. One key decision is the interface for creating graphics. At this point there are at least three different ones (Winston, Gadfly, and Gaston), and perhaps more will be generated before a dominant one is arrived at. As such, we don't try to teach the details of any one of them. That being said, we can do all we need with the easy-to-use plot function from the Gadfly package.
This function makes graphing a function as easy as specifying the domain to graph over, e.g.:
f(x) = exp(-x^2/2)
plot(f, -3, 3) Within each session you must first bring the package features into your workspace with the using command. Here is how we load the Gadfly package which will draw our graphics:
using GadflyIf you wanted to use Winston's graphics, then the command would be using Winston.
The most basic usage for plotting a function follows this pattern:
plot(function_object, from, to)as in
plot(sin, 0, 2pi)That creates the graphic. If you are working within IJulia it will be automatically displayed as an SVG file that allows you to zoom and pan with the mouse. The size of the graphics produced can be adjusted through a command like: set_default_plot_size(9inch, 6inch).
If working from a command line, Gadfly plots will be rendered in a browser. If desired, plots can also be written to a file using the draw function.
The above example used sin a function object. We can also use anonymous functions which are useful for one-off usage:
plot(x -> exp(x) + exp(-x), -2, 2) The last two examples show the two styles you can use to specify the function object. Either you have a name for the function object, such as sin, of 'f' (with no (x)!) or an anonymous function.
Again, we plot a function, this time a basic polynomial:
f(x) = x^2 - 2x + 2
plot(f, -3, 3) Using an anonymous function, we now graph that same function squared:
plot(x -> f(x)^2, -3, 3)Gadfly plots functions by creating about 250 paired points \((x, f(x))\), plots these, then connects the points with line segments. In the 250 function evaluations, it is of course quite possible that some of the points will return Inf or NaN.
The values which are Inf can not reasonably be plotted and when encountered Gadfly will issue a BoundsError, as in the value is out of bounds.
Values which are NaN are treated differently. Such points are simply not plotted, and no line segments are drawn causing the plot to be discontinuous. This convention can be utilized to good effect.
For example, this use truncates large values of a function:
trim_large(f; val=10) = x -> abs(f(x)) > val ? NaN : f(x)
f(x) = 1/x
plot(trim_large(f), -3, 3)This function trim_large returns a function that just replaces any large values (in absolute value) of f with NaN. This also works with Inf values, as Inf can be ordered by >.
Plot the function \(f(x) = x^3 - x\) over \([-2,2]\). How many zeros are there?
Plot the function \(f(x) = x^3 - x\). When is the function positive?
Plot the function \(f(x) = 3x^4 + 8x^3 - 18x^2\). Where (what \(x\) value) is the smallest value?
Plot the function \(f(x) = 3x^4 + 8x^3 - 18x^2\). What is the smallest value?
Plot the function \(f(x) = 3x^4 + 8x^3 - 18x^2\). When is the function increasing?
Plot the polynomial function \(p(x) = x^3 - 7.11x^2 + 16.8503x - 13.311105\). How many real roots are there \(1\) or \(3\)?
A rational function is nothing more than a ratio of polynomial functions, say \(f(x) = p(x)/q(x)\). One interesting thing about them is there can be asymptotes. These can be vertical (which can happen when \(q(x)=0\)), horizontal (as \(x\) gets large or small), or even slant.
The vertical asymptotes require care when plotting, as the naive style of plotting where a collection of points is connected by straight lines can render poor graphs when the scale of \(y\) values get too large. As well, the plot function does not handle cases where you get Inf for an answer (values that can come up when division by 0 is possible).
Keep in mind: while a graphic can highlight different features of a graph, it is often not possible to look at all of them on one graph.
Some features of a graph that are identifiable by calculus concepts are:
For example, if you want to find zeroes of a function, you really want to look at areas of the graph where the \(y\) values are close to \(0\). However, if you have a vertical asymptote on the same graph, the \(y\) values might also be asked to show very large or small values. With only a finite number of pixels available, it is impossible to easily do both.
What to do? If you were on a smartphone, you might be tempted to pan around to avoid the asymptotes, then pinch and zoom narrow the graph to the feature of interest. With julia you basically do the same thing, though panning and zooming is done by changing the domain of the \(x\) values used in plotting.
Let's look at the function \(f(x) = 1/x\), which has an obvious vertical asymptote at \(0\).
One can try a simple plot:
plot(x -> 1/x, -3, 3) The issue at \(0\) is avoided, as the points chosen by Gadfly do not include \(0\). The asymptote appears as a strongly slanted line, as individual points are simply connected in a dot-to-dot manner.
Doing better without much work is done by simply restricting the part that is plotted. (Alternatively, you can use NaN values or multiple functions on one.) For this example, we can easily draw the positive part:
It is best to avoid the asymptote at \(0\) completely by backing away by enough to avoid the large range on the \(y\) axis. In this case, starting around \(1/10\) is reasonable:
plot(x -> 1/x, 1/10, 3)Let's look at the rational function
\[ f(x) = \frac{(x-2)(x-3)}{x-1} \]
This will have a vertical asymptote at \(1\), zeros at \(2\) and \(3\) and a slant asymptote. A quick graph from \(-10\) to \(10\) shows just some of this:
f(x) = (x-2)*(x-3)/(x-1)
plot(f, -10, 10)We can see the slant asymptote and the vertical asymptote, but have no chance of seeing the zeros of the local extrema. We can if we restrict our plot range.
For example, we graph to the right of the asymptote with:
plot(f, 1 + 0.5, 4)This shows the two zeros and gives an idea of the relative minimum. Similarly, a plot of the left of the asymptote can be illustrative. Here we step back by a bit:
plot(f, -3, 1 - 0.1) It appears that the relative maximum occurs between \(-1\) and \(0\).
Here we see what happens to the asymptote. The scale of the \(y\) values is huge. We added a small amount to the left endpoint in case the function is not defined there, but this function takes the reciprocal of a small amount and makes it huge. Clearly we need to really avoid the issue. It isn't hard -- just add a little bit more to \(0\):
Example:Trigonometric functions
Let
\[ f(x) = \frac{5}{\cos x} + \frac{10}{\sin x}. \]
Estimate graphically the minimum value over \((0, \pi/2)\).
The domain comes from the fact that \(\sin(0) = 0\) and \(\cos(\pi/2) = 0\), so we will have asymptotes at each. A simple graph shows there are issues:
As typical with vertical asymptotes, we can't see finer features of the graph when the asymptotes are drawn. To remedy, we again back off from the boundaries. Since \(\sin(x)\) behaves like \(x\) near \(0\), we pick delta = 0.3 again and expect a max near \(10/(3/10) \approx 33\).
delta = 0.3;
plot(f, 0 + delta, pi/2 - delta) With this, we see the minimum value is near \(y=20\) and occurs around \(0.9\).
The function \(f(x) = (x^3 - 2x) / (2x^2 -10)\) is a rational function with issues when \(2x^2 = 10\), or \(x = -\sqrt{5}\) or \(\sqrt{5}\).
Plot this function from \(-5\) to \(5\). How many times does it cross the \(x\) axis?
The function \(f(x) = (x^3 - 2x) / (2x^2 -10)\) has a slant asymptote. One can find the slope of this line using algebra, but if you prefer the computer, you can graph. Define both
f(x) = (x^3 - 2x) / (2x^2 -10)
g(x) = abs(x) > 5 ? f(x) : NaNthen plot g over \([-20, 20]\). Using algebra or this graph, estimate the slope?
The rational function \(f(x) = (x^2 - 3x - 2) / (x^2 - 4)\) has issues at \(-2\) and \(2\). How many times does its graph cross the \(x\) axis?
When we learn how to make a graph using paper and pencil, the "T" method is employed, so called as we draw a "T" and fill in values to plot for \(x\) and \(y\).
For example, a chalkboard after the instructor shows how to plot \(f(x) = x^2\) might have this drawn on it:
x | y
------
1 | 1
2 | 4
3 | 9
4 | 16We would like to be able to mimic the following procedures used above in julia:
As these are the steps done to create the ordered pairs for a plot.
We have seen how variables can be used to refer to a single value, but we want to refer to more than one value now. A container for holding such is an Array. Arrays are implemented in most all computer languages, though the term can mean different things. We are looking for vectors, or a one-dimensional array. In general, an array is a multidimensional grid of values -- all of the same type (integer, floating point, functions, ..., or ANY).
For our purposes, we want vectors (one dimensional, \(n\) by 1 arrays in julia). These can be easily constructed in different ways.
[]One can directly create a 1-dimensional array with the [] syntax. For example, to put a sequence of numbers together we can do:
[1, 2 , 3, 4, 5]
## [1,2,3,4,5]Or
[1, 1, 2, 3, 5, 8]
## [1,1,2,3,5,8]These create "vectors." (Our printing in the notes displays these horizontally to save space, but the output would appear vertically if done at the console.) Row vectors (which are arrays, but not of type Vector) can be created by dropping the commas:
[13 21 34 55]
## 1x4 Array{Int64,2}:
## 13 21 34 55The notation is that [a,b,c] combines a, b, and c vertically (vcat) and [a b c] combines them horizontally (hcat). The former is good to make column vectors for single values (scalars).
Row and column vectors are different! We will use just column vectors going forward.
A basic set of numbers used in programming are the values 1,2, ..., n. These are the simplest example of an arithmetic progression or sequence, which in general can start in a different place and have steps of a size different from \(1\):
\[ a, a + h, a+2h, a+3h, ..., a + nh \]
It should be possible to specify arithmetic sequences either by
the start and end points and the number of points employed
the start and end points and step size.
In julia the linspace function will do the former and the range operator the latter.
Here are 5 evenly spaced numbers from \(0\) to \(\pi/2\) given by linspace (linearly spaced numbers):
linspace(0, pi/2, 5)
## [0.0,0.39269908169872414,0.7853981633974483,1.1780972450961724,1.5707963267948966]A column vector is returned.
The 5 above is optional -- leaving it out will specify the default of \(100\) points:
linspace(0, pi/2)
## [0.0,0.015866629563584814,0.03173325912716963,0.04759988869075444,0.06346651825433926,0.07933314781792407,0.09519977738150888,0.11106640694509369,0.12693303650867852,0.14279966607226333 … 1.427996660722633,1.443863290286218,1.4597299198498028,1.4755965494133878,1.4914631789769726,1.5073298085405573,1.523196438104142,1.5390630676677268,1.5549296972313118,1.5707963267948966]The "range" operator, :, is used to specify a step size, like \(h\) in the definition above. To get values which step by 1 we simply specify the start and end values:
1:4
## 1:4That isn't so impressive. The description julia uses to show this value is exactly how we defined it, but this range is specifying the values 1, 2, 3, 4. To see that, we can put inside, [] to make an array:
[1:4]
## [1,2,3,4]The range operator returns a Range object which internally is much more compact to store than the array, though at times less useful. (At those times, they can be converted to arrays as above.)
The range operator can also work with a step size:
a=0; b=10; h=3
a:h:b
## 0:3:9
[a:h:b]
## [0,3,6,9]Notice, the value for b is treated as as suggestion, the range will stop without exceeding b.
The \(x\) values for a plot are typically a sequence of increasing values from \(a\) to \(b\). We would generally like to be able to specify the number of values to plot. This makes linspace the go-to choice to use. (Though Gadfly uses: a:(b-a)/n:b with n=250.)
Which command will produce the sequence \(1,3,5,7,9,11, ..., 99\)?
Which command produces 10 numbers between 0 and 10 that are evenly spaced?
Which command does not produce the numbers \(1, 2, 3, 4, 5\)?
Which command does produces the numbers \(1, 1, 2, 3, 5, 8, 13\)?
A column vector has a natural sense of first, second, ..., the \(n\)-th element. This allows julia to refer to the values by index (\(1\)-based, unlike some other computer languages). So, if x is an array, x[1] is the first value in that array. One can extract and assign values using indices. A simple example is:
x = [2, 4, 6, 8]
## [2,4,6,8]
x[1]
## 2
x[3]
## 6There are some special values. The end value refers to the last (\(n\) th):
x[end]
## 8The \(n\) -- or number of elements -- can be returned by length:
length(x)
## 4A range object can be used for indices:
x[1:3]
## [2,4,6]As can any other column vector (but not a row vector):
x[ [ 1, 2, 3] ]
## [2,4,6]The value end can be used in a range when indexing:
x[2:end]
## [4,6,8](But not without indexing, as you can see by typing 2:end by itself.)
Let
x = [1, 1, 2, 3, 5, 8, 13]
## [1,1,2,3,5,8,13]What is the value of x[end - 1] + x[end]?
Let
x = [1, 1, 2, 3, 5, 8, 13]
## [1,1,2,3,5,8,13]What is the value of x[3]?
When a vector is created, if possible the resulting values are converted to be the same type. Let
x = [1, 2.0]
## [1.0,2.0]
y = [1, 2.0, "three"]
## {1,2.0,"three"}For x[1] and y[1] what does typeof return?
(The y container is of type Any which allows it to hold any type of object, the x container only holds values of a certain type.)
To specify the \(y\) values we wish to "map" the function f to each \(x\) value. In julia there are many different ways to do this, we list four.
map function. In many areas of mathematics, one refers to a function as a "map" from the domain to the range. The implication is that the function takes all the \(x\) values to the corresponding \(y\) values at once (conceptually) and not one at a time. The map function will apply the function f to each value in the array x, basically taking [x1, x2, ..., xn] and returning [f(x1), f(x2), ..., f(xn)].For example, let's look at the simple polynomial \(f(x) = -16x^2 + 32x\). We define our julia function with:
f(x) = -16x^2 + 32xIf we want to look at this function for \(x\) values between \(0\) and \(2\) we might define the \(x\) values with:
x = linspace(0, 2, 5)
## [0.0,0.5,1.0,1.5,2.0]Then the map function will create the corresponding \(y\) values:
map(f, x)
## [0.0,12.0,16.0,12.0,0.0]The syntax of map requires a slight pause. Here we do not actually call the function f, as in f(2). Rather, we pass the name of the function object to the map argument -- and map calls the function for each value in the column vector x and returns a corresponding column vector.
It is also quite common to use anonymous functions with map. For example:
map(u -> 10 * u^2, x)
## [0.0,2.5,10.0,22.5,40.0]We use u for the dummy variable in the anonymous function, so as not to get it confused with the variable x holding our values, but this is not necessary.
\[ y = \{ f(x): x \text{ in } [0, 2] \}, \]
where this is read the values \(f(x)\) for each \(x\) in the interval \([0,2]\).
Here we define values vals to represent the continuum of values in the interval \([0,2]\), then use a "comprehension" to create the set notation above. The syntax is similar:
vals = linspace(0, 2.0, 5)
## [0.0,0.5,1.0,1.5,2.0]
[f(x) for x in vals] ## or {f(x) for x = vals}
## {0.0,12.0,16.0,12.0,0.0}Square brackets are typical when using a comprehension, though curly braces are possible -- try it. While curly braces are more mathematically familiar -- they always return a container of type Any. Square brackets try to guess a better type, if possible.
In this example we still get a container of type Any with square brackets. While not critical, we would really like to have this be of type "float", as some math functions care about this distinction. Converting can be done a few ways:
float (which coerces the values):float( [ f(x) for x in vals ] )
## [0.0,12.0,16.0,12.0,0.0]Float64[f(x) for x in vals]
## [0.0,12.0,16.0,12.0,0.0]Either way, perhaps map is a bit less trouble. (For ranges it is just syntactic sugar, as a map is defined by a comprehension.) Regardless, we mostly use comprehensions going forward as they mirror a more familiar mathematical syntax and generalize to functions of more than one variable nicely.
for loop: The for loop is a very common pattern in computer programming. For speed reasons, some languages (e.g., MATLAB and R) try to avoid for loops for a "vectorized" approach (see below), but julia works well with for loops, and they are generally faster than a vectorized approach.A for loop simply loops over each value in the iterable data vector x giving it a temporary name as it goes. To use a for loop to create the \(y\) values requires several steps. The first creates a container to hold the new values, using the similar function to make a vector of the same size as x below (but fills it with nonsense):
y = similar(x);
for i in 1:length(x)
y[i] = f( x[i] )
endThe for loop ends with the keyword end. Here we loop over each index of x and assign to the corresponding y value f(x[i]).
Conceptually this is the opposite of map where we think of the function acting on the entire column vector x. Instead, we iterate one-by-one over the values of x saving the function values as we go. The use of for loops is called imperative programming, as you describe each step. The use of functions like map is called declarative programming as you simply declare what you want and the computer language does the rest.
In some languages for loops are avoided if possible, as they can be slower. As well, they can require extra bookkeeping, such as needing to allocate a container for the answers. That being said, in julia they are widely used for speed and storage reasons. As well, they are used when we need to refer to more than one index. An example of that is the following way to create a Fibonacci pattern from the formula \(x_i = x_{i-1} + x_{i-2}\):
x = zeros(Int, 10); ## pre-allocate an integer array
x[1:2] = [1,1]
## [1,1]
for i in 3:length(x)
x[i] = x[i-1] + x[i-2]
endRelated to a for loop is the while loop. This will repeat as long as some condition is true. The following pattern reproduces a for loop:
i, n = 1, length(x)
## (1,10)
while (i <= n)
print( x[i], " " ) ## do something ...
i = i + 1
endWe will use while loops when we iterate some process and are waiting until some computed value gets close to \(0\).
julia allows for "vectorization" of values. That means, many functions work on an array all at once. Manysin(x)
## [0.8414709848078965,0.8414709848078965,0.9092974268256817,0.1411200080598672,-0.9589242746631385,0.9893582466233818,0.4201670368266409,0.8366556385360561,0.5290826861200238,-0.9997551733586199]but not all:
x^2
## "MethodError(*,([1,1,2,3,5,8,13,21,34,55],[1,1,2,3,5,8,13,21,34,55]))"The latter fails, whereas the sin function naturally returns the 5 values for x. The reason is the definitions of multiplication, division and powers are mathematically different for single numbers (scalars) than for column vectors, or more generally matrices. Like MATLAB, julia chose the simplest notation to do matrix multiplication, powers, and division. To get element-by-element operations, we prefix the operator with a "dot", as in .^:
x .^ 2
## [1,1,4,9,25,64,169,441,1156,3025]So, if we were to redefine our function \(f(x) = -16x^2 + 32x\) with:
f(x) = -16x.^2 + 32xThen we could have simply done:
f(x)
## [16,16,0,-48,-240,-768,-2288,-6384,-17408,-46640]Though this is easy to do, it is also rather tedious, as putting in the "dots" is, for some reason, not that intuitive, so mistakes are often made. (It can surprisingly, also be slower to perform.)
We do use this style to generate simple values, for example here we step through some powers of \(1/10\) without much fuss:
(1/10).^(1:5)
## [0.1,0.010000000000000002,0.0010000000000000002,0.00010000000000000005,1.0000000000000003e-5]Putting this altogether, to create the "T"-table used to graph \(y=x^2\), we could do any of these:
begin
x = [1:4]
y = map(u -> u^2, x)
[x y] ## integer type
end
## 4x2 Array{Int64,2}:
## 1 1
## 2 4
## 3 9
## 4 16or
begin
x = [1, 2, 3, 4]
y = [x^2 for x in x]
[x y] ## Any type, as y is -- not the same
end
## 4x2 Array{Any,2}:
## 1 1
## 2 4
## 3 9
## 4 16or
begin
x = 1:4
y = similar(x)
for i in 1:length(x)
y[i] = x[i]^2 ## integer, as x[i] is
end
[x y]
end
## 4x2 Array{Int64,2}:
## 1 1
## 2 4
## 3 9
## 4 16begin
x = linspace(1.0, 4.0, 4)
y = x.^2
[x y] ## Float type, as x is here
end
## 4x2 Array{Float64,2}:
## 1.0 1.0
## 2.0 4.0
## 3.0 9.0
## 4.0 16.0The comments show that each is slightly different for technical reasons, but all display the same values.
Does this command produce the values \(\{.1, .01, .001, .0001\}\)?
x = [(1/10)^i for i in 1:4];Does this command produce the numbers \(10, 100, 1000, 10,000\)?
map(i -> 10^i, 1:4);Let \(f(x) = x^2 - 2x\). Which command will produce a \(y\) values for plotting \(f\) over \([0, 2]\)?
Will this command produce \(y\) values for plotting \(f(x)\) over \([0,1]\)?
[f(x) for x in 0:1/100:1];Sometimes the function is used to ask a question. For example, this command answers if the values are prime or not:
[isprime(i) for i in 1:10]
## [false,true,true,false,true,false,true,false,false,false]This shows there are 4 primes in 1 to 10. How many are there between 100 and 110?
(It is easier to write map(isprime, 1:10) in this case.)
If one has two vectors, x and y of equal size, then creating a graphic for them is straightforward. The basic syntax is
plot(x=x, y=y, Geom.line)For example, to plot \(y=x^2\) over \([-1,1]\) we might do:
f(x) = x^2
xs = linspace(-1, 1)
ys = map(f, xs)
plot(x=xs, y=ys, Geom.line)The specification Geom.line instructs Gadfly how to present the data, in this case connected with lines. One can add Geom.point to get the points filled in as well. This is the default setting, so if omitted a scatter plot is drawn.
Now that we know about vectors, we can use them to illustrate how to graph two (or more) functions at once over the same interval -- just make a vector of functions and hand this off to plot.
For example, to graph both the sine and cosine function we have:
plot([sin, cos], 0, 2pi) Or to compare the effects of a simple transformation:
f(x) = x^2
g(x) = 4 + f(x-3)
plot([f, g], -4, 4)To print a heavier \(x\)-axis, we can graph the function \(y=0\), specified through the anonymous function x -> 0:
f(x) = x^2 - 2
plot([f, x -> 0], -2, 2)Similarly, this example adds a straight line at the origin, but only when the cosine is positive:
plot([sin, x -> cos(x) > 0 ? 0.0 : NaN], 0, 2pi) ## need 0.0 and not simply 0 hereNot all pairs of functions will naturally appear on the same \(x\) and \(y\) axes, as the scale of them may be way off. In fact that is typical. Yet, you may still wish to compare them. The vstack command in the Compose package (this package does the drawing for Gadfly) allows us to vertically stack graphics. (There is also hstack for horizontal stacks.)
The simplest usage just has one add graphics with commas:
using Compose
begin
p1 = plot(x -> x^2 + 3x + 10, -4, 4);
p2 = plot(x -> 2x + 3 , -4, 4);
vstack(p1, p2)
endThe vstack function stacks plots vertically. There is an accompanying hstack command for putting graphics side by side.
Define \(f(x)\) to be a triangular function as follows:
f(x) = abs(x) > 1 ? 0.0 : 1.0 - abs(x)In many applications, the following transformation is employed:
\[ g(x, c, h) = \frac{1}{h} f(\frac{x - c}{h}) \]
For constants \(h\) and \(c\).
Make a graph of both \(f(x)\) and \(g(x, 1, 1/2)\) over the interval \([-2,3]\). Consult the graph to see which statement is true?
We saw that this command will produce two graphs:
plot([sin, x -> cos(x) > 0 ? 0 : NaN], 0, 2pi)What is the sine curve doing when the flat line is drawn?
Make a graph of \(f(x) = x\), \(g(x) = \tan(x)\), and \(h(x) = \sin(x)\). Over the interval \([0,\pi/4]\). Based on this graph which of the following below seems correct?