
This page will walk you through basic numerical programming in julia.
Table Of Contents 

1. Running julia 
2. Simple programs on the fly 
3. More complicated programs: Program files 
4. Plotting in julia 
5. Plotting in gnuplot 
6. Arrays, vectors and matrices 
One way to use julia is to work at an Athena cluster machine. To get it running, use the following command:
username@athena: ~$ athrun julia julia
Unfortunately, the version of julia
installed on most of the Athena dialup machines
(such as linux.mit.edu
) is out of date and
broken. The exception is test.dialup.mit.edu,
on which the above command (athrun julia julia
)
works. However, in any event you should not be using
dialup machines for computationintensive programs; use
Athena cluster machines instead.
Alternatively, it is quite easy to install julia on your own computer. (Indeed, one of the key virtues of the language is that it is free!) However, this page does not cover installation; for instructions on how to install julia on your Linux, Windows, or Mac system, see here or here. For personalized handson assistance, come to office hours or email me to arrange a time to meet.
For now, we assume you have installed and run julia, so that you should see something like this on your screen:
For quick and easy calculations, you can type commands right into the julia command prompt.
julia> tanh(3.4) + sinh(5.6) 136.2091297091523 julia> exp(pi) 23.140692632779267
Here's an approximate evaluation of the sum we considered on the first day of class.
julia> Sum=0.0 julia> for n=1:2:10000 Sum += 1/n + 1/(n+1); end julia> Sum 0.6930971830599458 julia>
function PartialSum(N) Sum=0.0 for n=1:2:N Sum += 1/n + 1/(n+1) end Sum end
Save this to a file named PartialSum.j
and load it into your julia
session using the require
command:
julia> require("PartialSum.j")
Now we can repeat the calculation we did onthefly above like this:
julia> PartialSum(10000) 0.6930971830599458
To investigate the convergence of the partial sums, let's
evaluate our PartialSum
function at 20
values of N between 10^{2} and 10^{8}:
julia> NVector=int(logspace(2.0,8.0,20)) 20element Int64 Array: 100 207 428 ... 23357215 48329302 100000000
The command logspace(a,b,n)
creates
a logarithmicallyspaced vector of n
numbers between
10^{a}
and
10^{b}.
Then we use int
to round each element of
this vector to the nearest integer.
julia> SVector = 0.0*NVector; julia> for n=1:length(NVector) SVector[n] = PartialSum( NVector[n] ); end
Here we first defined SVector
to be a
vector of the same length as NVector
,
then filled in each element of SVector
with the value of PartialSum
evaluated
at the corresponding element of NVector.
To get a sidebyside visualization of the results, we can go like this:
julia> [ NVector SVector ] 20x2 Float64 Array: 100.0 0.688172 207.0 0.690749 428.0 0.69198 886.0 0.692583 1833.0 0.692875 3793.0 0.693015 7848.0 0.693083 16238.0 0.693116 33598.0 0.693132 ... 297635.0 0.693146 615848.0 0.693146 1.27428e6 0.693147 2.63665e6 0.693147 5.4556e6 0.693147 1.12884e7 0.693147 2.33572e7 0.693147 4.83293e7 0.693147 1.0e8 0.693147 julia>
The command [NVector SVector]
concatenates
our two 20element vectors into a 20×2 matrix, which
is then printed to the console.
Based on this, it looks like the 6th digit of the partial sum is stabilizing somewhere around N=10^{6}, but let's be more precise.
Recall that the exact (N→∞) value of our sum
is ln(2). Let's obtain a vector of the same size as
SVector
whose nth element gives us
the error incurred by truncating the partial sum at
N=NVector[n]
.
julia> ExactVector=log(2.0) * ones( length(NVector) ) 20element Float64 Array: 0.693147 0.693147 0.693147 ... 0.693147 0.693147 0.693147 julia> AbsErrorVector=abs(ExactVector  SVector) 20element Float64 Array: 0.004975 0.00239807 0.00116686 ... 2.14035e8 1.03558e8 5.08625e9
In the first line here, we used the ones
command to obtain a vector of the same length as
NVector
in which every entry is 1. Then
we multiplied that vector by log(2.0)
to
obtain a vector whose every entry is equal to the exact
value of our sum. Finally, we defined AbsErrorVector
to be a vector whose nth entry is the absolute
error between the exact answer and the partial sum
evaluated at N=NVector[n].
It's also convenient to look at the relative error, which is the absolute error divided by the magnitude of the exact answer:
julia> RelErrorVector=AbsErrorVector / log(2.0) 20element Float64 Array: 0.00717741 0.00345968 0.00168342 ... 3.08788e8 1.49402e8 7.3379e9
In this case the exact answer is on the order of 1, so the relative and absolute errors don't differ much. But they are quite different quantities when the exact answer we are computing is a large or a small number.
As of February 2013, plotting is not part of the core julia system, but instead is handled by a separate addon package named winston. Before you can follow the examples discussed below, you will need to make sure the winston package is present within your julia installation. Instructions for installing winston are discussed here. (There is also some incomplete documentation at the winston github page.)
Once you have successfully added the winston package to your julia installation, you will want to run the following command each time you start a julia session:
julia> using Winston
This loads all the winston commands into your session so that they are available for you to run at the julia command line.
julia, like matlab, doesn't really plot functions  it plots vectors. To create a plot, you need to create two vectors, with equal lengths; one vector will contain the xcoordinates of the data points plotted, and the other vector contains the y coordinates. For example, to plot sin(x) between 0 and 2π, you can go like this:
julia> xVector=[0:0.01:2*pi]; julia> yVector=0.0*xVector; julia> for n=1:length(xVector) yVector[n] = sin(xVector[n]); end julia> plot(xVector, yVector)
For more complex plotting tasks such as adding axis
labels, plotting multiple curves on a plot, and
saving your plot to a file, you will create a Julia
plot object using the FramePlot()
function. This routine returns a bare plot object,
to which you add curves and annotations as you like.
When you are finished constructing your plot object,
you use the file
command to save it to
a file.
xVector=[0:0.01:2*pi] sinVector=sin(xVector); cosVector=cos(xVector); sinCurve=Curve(xVector, sinVector); setattr(sinCurve, "label", "Sine curve"); cosCurve=Curve(xVector, cosVector); setattr(cosCurve, "label", "Cosine curve"); MyLegend=Legend(0.5, 0.9, {sinCurve, cosCurve} ) p=FramedPlot() setattr(p,"title","Sine and cosine curves"); setattr(p,"xlabel","X data"); setattr(p,"ylabel","Y data"); add(p, sinCurve, cosCurve, MyLegend); file(p, "SineCosine.pdf"); file(p, "SineCosine.png");
To generate loglinear or loglog plots, you can assign one or more of the axes to be logarithmic like this:
setattr( p.x1, "log", true ) # make lower xaxis logarithmic setattr( p.x2, "log", true ) # make upper xaxis logarithmic setattr( p.y1, "log", true ) # make left yaxis logarithmic setattr( p.y2, "log", true ) # make right yaxis logarithmic
You can find some more examples of winston plotting here.
gnuplot is an outstanding free opensource software tool that you may use to produce highquality plots. The simplest way to plot a data curve in gnuplot is to write the data you want to plot into a text file, with the numbers for the xaxis and yaxis data in separate columns separated by spaces. Thus, using gnuplot to do plotting in julia is a twostep process:
In this case, we want to plot the entries of NVector
on the xaxis and the entries of AbsErrorVector
or RelErrorVector
on the yaxis. To write these
data to a file of the correct format, we'll use
julia's
writedlm
function, which writes a matrix to a text file.
This function takes three arguments,
writedlm(FileName, Matrix, Separator)

where the meaning of the arguments is:
FileName
: a string like "Error.dat"
(note the double quotes) that specifies the name of the file
we want to write.
Matrix
: a julia
matrix. If the matrix has dimension M×N,
then the output file will have M lines, each
containing N numbers.
In our case, we actually have three vectors of length
20, NVector
, AbsErrorVector
,
and RelErrorVector
.
We can concatenate these into a matrix of size
20×3 by saying
NewMatrix=[NVector AbsErrorVector RelErrorVector]
In the code below, we will pass the righthandside of this
construction directly as the Matrix
argument to writedlm.
Separator
: a character like
';'
or ','
(note the single quotes) that will be used to separate
the columns of numbers from each other on each line
of the output file. Here we will choose this to be
simply a blank space, ' '
.
After executing the following julia
code:
writedlm( "Error.dat", [NVector AbsErrorVector RelErrorVector], ' ')
we have a file named Error.dat
with the
following content:
100 0.004975001249750366 .0071774096314312485 207 0.002398067744293408 .003459680442407883 428 0.001166859554777133 .0016834224930908739 ... 23357215 2.14035407175572e08 3.08787820506848e08 48329302 1.03557514785635e08 1.49401913028021e08 100000000 5.08624675710223e09 1.03557514785634e08
The first column is N, the second column is the absolute error in the partial sum S_{N}, and the third column is the relative error in S_{N}.
In this case, the x axis data we want to plot
are in column 1 of Error.dat
, and the yaxis
data are in columns 2 and 3. To plot column 2 vs. column 1,
launch
gnuplot and type the following
commands:
556 maki /home/homer <> gnuplot gnuplot> set logscale xy gnuplot> plot 'Error.dat' u 1:2 gnuplot>
The command set logscale xy
requests a loglog
plot. The command plot 'Error.dat' u 1:2
says "plot the contents of file Error.dat
with
column 1 as the x axis and
column 2 as the y axis." (The u
stands
for "using.") This should produce something like the
following:
To make this look slightly nicer, we can add a plot of the relative error, add labels for the axes and a title for the graph, and connect the data points with lines:
gnuplot> set xlabel 'N' gnuplot> set ylabel 'Error in partial sum S_N' gnuplot> set title 'Convergence of partial sums' gnuplot> set logscale xy gnuplot> plot 'Error.dat' u 1:2 t 'Abs Error' w lp, '' u 1:3 t 'Rel Error' w lp
This produces something like this:
(The option w lp
to the plot
command stands for "with lines and points". This
instructs
gnuplot
to plot both data points and a line through the
points. You can experiment with various other
gnuplot options
to customize your plots to your taste.)
Once you have configured all your gnuplot options just right to obtain a beautiful plot on your screen, you can capture this plot to a postscript file for printing or for inclusion in other documents using the following commands:
gnuplot> set terminal postscript color gnuplot> set output 'MyPlot.eps' gnuplot> replot
This produces a file named MyPlot.eps
containing a postscript rendition of the last plot
you generated with the 'plot' command. You may
then open this file in a postscript viewer, print it,
convert it to PDF, include it in LaTeX documents, etc.
Each of the following commands creates a vector of length 10:
julia> xVector = [1:10] julia> hyVector = ones(10); julia> zVector = 7.8 * yVector;
After creating a vector, you can access its nth element using square brackets (not parentheses as in matlab):
julia> q = zVector[3]; # sets q = 7.8 julia> zVector[3] = 4.5; # changes the value of zVector[3]
Each of the following commands creates a matrix of dimension 10×20:
julia> xMatrix = zeros(10,20); julia> yMatrix = 2.3 * xMatrix;
As before, you use square brackets to access the components of matrices:
julia> q = xMatrix[3,8]; # sets q = 0 julia> yMatrix[3,8] = 5.75; # changes the 3,8 entry of yMatrix
If you have two N dimensional vectors, you can concatenate them to form an N×2 matrix. To do this, just put the vectors sidebyside inside square brackets, as though you were aligning the columns of your matrix sidebyside:
julia> xVector = [1:10]; julia> yVector = [101:110]; julia> xyMatrix = [xVector yVector]; # creates a 10x2 matrix julia> xyMatrix[4,2] 104
You can also do the opposite of concatenation: you
can extract individual columns or rows of
matrices. The tool you want to use here is the colon
character (:
), which you can use in place
of a number in square brackets to mean "the whole row"
or "the whole column." For example, to unpack
xyMatrix
into xVector
and
yVector
in the above code snippet we
could say
julia> xVector = xyMatrix[:,1]; # set xvector = entire 1st column of xyMatrix julia> yVector = xyMatrix[:,2]; # set yvector = entire 2nd column of xyMatrix
It's also possible to define vectors obtained as individual rows of a matrix. For this purpose we put the colon character as the second index in the square brackets (meaning "give me the entire row"):
julia> ThirdRow = xyMatrix[3,:] # get entire 3rd row of xyMatrix 1x2 Int64 Array: 3 103 julia>
18.330 Spring 2013  Basic Numerical Programming in julia 