[an error occurred while processing this directive]

MIT Home Page

18.330 Introduction to Numerical Analysis

Basic Numerical Programming in Julia

Spring 2013

Instructor: Homer Reid

MIT Home Page


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

Running julia

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 computation-intensive 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 hands-on 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:

Simple programs on-the-fly

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>

More complicated programs: Program files

For more complicated situations, you will probably want to define your own julia functions in files, which you can save and re-use in multiple julia sessions. For example, to create a function that evaluates the partial sum above with a user-specified value of N, open your favorite text editor and enter the following program:

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 on-the-fly 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 102 and 108:

julia> NVector=int(logspace(2.0,8.0,20))

20-element Int64 Array:
       100
       207
       428
       ...
  23357215
  48329302
 100000000

The command logspace(a,b,n) creates a logarithmically-spaced vector of n numbers between 10a and 10b. 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 side-by-side 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 20-element 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=106, but let's be more precise.

Looking at the error vector

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) )
20-element Float64 Array:
 -0.693147
 -0.693147
 -0.693147
  ...
 -0.693147
 -0.693147
 -0.693147

julia> AbsErrorVector=abs(ExactVector - SVector)

20-element Float64 Array:
  0.004975   
  0.00239807 
  0.00116686 
  ...
  2.14035e-8 
  1.03558e-8 
  5.08625e-9 

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)
20-element Float64 Array:
 0.00717741 
 0.00345968 
 0.00168342 
 ...
 3.08788e-8 
 1.49402e-8 
 7.3379e-9  

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.

Plotting in julia

As of February 2013, plotting is not part of the core julia system, but instead is handled by a separate add-on 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.

Simple plotting

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 x-coordinates 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)

More complex plots

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 log-linear or log-log plots, you can assign one or more of the axes to be logarithmic like this:

 setattr( p.x1, "log", true ) # make lower x-axis logarithmic
 setattr( p.x2, "log", true ) # make upper x-axis logarithmic
 setattr( p.y1, "log", true ) # make left y-axis logarithmic
 setattr( p.y2, "log", true ) # make right y-axis logarithmic

You can find some more examples of winston plotting here.

Plotting using gnuplot

gnuplot is an outstanding free open-source software tool that you may use to produce high-quality 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 x-axis and y-axis data in separate columns separated by spaces. Thus, using gnuplot to do plotting in julia is a two-step process:

  1. Write the data we want to plot from julia to a text file.
  2. Plot the contents of that file in gnuplot.

1. Writing julia data to a text file

In this case, we want to plot the entries of NVector on the x-axis and the entries of AbsErrorVector or RelErrorVector on the y-axis. 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:

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.14035407175572e-08  3.08787820506848e-08
 48329302  1.03557514785635e-08  1.49401913028021e-08
100000000  5.08624675710223e-09  1.03557514785634e-08

The first column is N, the second column is the absolute error in the partial sum SN, and the third column is the relative error in SN.

2. Plotting a text file in gnuplot:

In this case, the x- axis data we want to plot are in column 1 of Error.dat, and the y-axis 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 log-log 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.)

Creating postscript output from gnuplot plots

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.

Arrays, vectors, and matrices

julia shares with matlab the ability to manipulate vectors and matrices conveniently.

Creating vectors

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]

Creating matrices

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

Concatenating: assembling vectors into matrices

If you have two N dimensional vectors, you can concatenate them to form an N×2 matrix. To do this, just put the vectors side-by-side inside square brackets, as though you were aligning the columns of your matrix side-by-side:

  julia> xVector = [1:10];
  julia> yVector = [101:110];
  julia> xyMatrix = [xVector yVector];   # creates a 10x2 matrix
  julia> xyMatrix[4,2] 
    104

Indexing: extracting vectors from matrices

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