Introduction to GNU Octave

Hubert Selhofer, revised by Marcel Oliver
updated to current Octave version by Thomas L. Scofield


Date: 2008/08/16


\includegraphics[height=10cm]{sombr}


Contents

1 Basics

1.1 What is Octave?

Octave is an interactive programming language specifically suited for vectorizable numerical calculations. It provides a high level interface to many standard libraries of numerical mathematics, e.g. LAPACK or BLAS.

The syntax of Octave resembles that of Matlab. An Octave program usually runs unmodified on Matlab. Matlab, being commercial software, has a larger function set, and so the reverse does not always work, especially when the program makes use of specialized add-on toolboxes for Matlab.

1.2 Help!

Example:

  octave:1> help eig

1.3 Input conventions

1.4 Variables and standard operations

When debugging user-defined objects, the following commands are useful:

Examples:

  octave:1> x12 = 1/8, long_name = 'A String'
  x12 = 0.12500
  long_name = A String
  octave:2> sqrt(-1)-i
  ans = 0
  octave:3> x = sqrt(2); sin(x)/x
  ans = 0.69846
And here is a script doless, saved in a file named doless.m:
  one = 1;
  two = 2;
  three = one + two;
Calling the script:
  octave:1> doless
  octave:2> whos
  *** local user variables:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   real scalar           1      1  three
   wd   real scalar           1      1  one
   wd   real scalar           1      1  two

2 Vector and matrix operations

Matrices and vectors are the most important building blocks for programming in Octave.


2.1 Vectors

Examples:

  octave:1> x = 3:6
  x =
    3  4  5  6
  octave:2> y = 0:.15:.7
  y =
  0.00000  0.15000  0.30000  0.45000  0.60000
  octave:3> z = pi:-pi/4:0
  z =
  3.14159  2.35619  1.57080  0.78540  0.00000

2.2 Matrices

A matrix $ \displaystyle A = \begin{pmatrix}1 & 2  3 & 4
\end{pmatrix}$ is generated as follows.

  octave:1> A = [ 1 2; 3 4]
  A =
    1  2
    3  4
Matrices can assembled from submatrices:
  octave:2> b = [5; 6];
  octave:3> M = [A b]
  M =

    1  2  5
    3  4  6

There are functions to create frequently used $ m \times n$ matrices. If $ m=n$ , only one argument is necessary.

2.3 Basic matrix arithmetic

Examples:

  octave:1> A = [1 2; 3 4]; B = 2*ones(2,2);
  octave:2> A+B, A-B, A*B
  ans =
    3  4
    5  6
  ans =
    -1   0
     1   2
  ans =
     6   6
    14  14


2.4 Element-wise operations

While * refers to the usual matrix multiplication, .* denotes element-wise multiplication. Similarly, ./ and .^ denote the element-wise division and power operators.

Examples:

  octave:1> A = [1 2; 3 4]; A.^2   % Element-wise power
  ans =

     1   4
     9  16

  octave:2> A^2   % Proper matrix power: A^2 = A*A
  ans =

     7  10
    15  22

2.5 Indexing and slicing

Querying dimensions: Reshaping:

Examples:

  octave:1> A = [1 2 3; 4 5 6]; v = [7; 8];
  octave:2> A(2,3) = v(2)
  A =

    1  2  3
    4  5  8

  octave:3> A(:,2) = v
  A =

    1  7  3
    4  8  8

  octave:4> A(1,1:2) = v'
  A =

    7  8  3
    4  8  8

2.6 Solving linear systems of equations

2.7 Inverses, decompositions, eigenvalues

A lot of these commands support further options. They can be listed by typing help funcname.

2.8 Testing for zero elements

3 Control structures


3.1 Functions

Traditionally, functions are also stored in plain text files with suffix .m. In contrast to scripts, functions can be called with arguments, and all variables used within the function are local--they do not influence variables defined previously.

Example:

A function f, saved in the file named f.m.

  function y = f (x)
    y = cos(x/2)+x;
  end

Remark:

In Octave, several functions can be defined in a single script file. Matlab on the other hand, strictly enforces one function per .m file, where the name of the function must match the name of the file. If compatibility with Matlab is important, this restriction should also be applied to programs written in Octave.

Example with two function values:

A function dolittle, which is saved in the file named dolittle.m.

  function [out1,out2] = dolittle (x)
    out1 = x^2;
    out2 = out1*x;
  end
Calling the function:
  octave:1> [x1,x2]=dolittle(2)
  x1 = 4
  x2 = 8
  octave:2> whos
  *** currently compiled functions:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   user function         -      -  dolittle
  *** local user variables:
  prot  type               rows   cols  name
  ====  ====               ====   ====  ====
   wd   real scalar           1      1  x1
   wd   real scalar           1      1  x2

Obviously, the variables out1 and out2 were local to dolittle. Previously defined variables out1 or out2 would not have been affected by calling dolittle.

3.2 Global variables

global name declares name as a global variable.

Example:

A function foo in the file named foo.m:

  global N   % makes N a global variable; may be set in main file

  function out = foo(arg1,arg2)
  global N   % makes local N refer to the global N
  <Computation>
  end
If you change N within the function, it changes in the value of N everywhere.

3.3 Loops

The syntax of for- and while-loops is immediate from the following examples:

  for n = 1:10
    [x(n),y(n)]=dolittle(n);
  end

  while t<T
    t = t+h;
  end
For-loop backward:
  for n = 10:-1:1
  ...

3.4 Branching

Conditional branching works as follows.

  if x==0
    error('x is 0!');
  else
    y = 1/x;
  end

  switch pnorm
    case 1
      sum(abs(v))
    case inf
      max(abs(v))
    otherwise
      sqrt(v'*v)
  end

3.5 Functions of functions

Example:

Approximate an integral by the midpoint rule:

$\displaystyle \int_a^b f(x)   dx \approx \frac{b-a}{N} \sum_{j=0}^{N-1}
f \left(a + (j+\frac{1}{2}) \frac{b-a}{N}\right)
$

We define two functions, gauss.m and mpr.m, as follows:
  function y = gauss(x)
    y = exp(-x.^2/2);
  end

  function S = mpr(fun,a,b,N)
    h = (b-a)/N;
    S = h*sum(feval(fun,[a+h/2:h:b]));
  end
Now the function gauss can be integrated by calling:
octave:1> mpr('gauss',0,5,500)

3.6 Efficiency considerations

Loops and function calls, especially through feval, have a very high computational overhead. Therefore, if possible, vectorize all operations.

Example:

We are programming the midpoint rule from the previous section with a for-loop (file name is mpr_long.m):

  function S = mpr_long(fun,a,b,N)
    h = (b-a)/N; S = 0;
    for k = 0:(N-1),
      S = S + feval(fun,a+h*(k+1/2));
    end
    S = h*S;
  end
We verify that mpr and mpr_long yield the same answer, and compare the evaluation times.
  octave:1> t = cputime;
  > Int1=mpr('gauss',0,5,500); t1=cputime-t;
  octave:2> t = cputime;
  > Int2=mpr_long('gauss',0,5,500); t2=cputime-t;
  octave:3> Int1-Int2, t2/t1
  ans = 0
  ans = 45.250

3.7 Input and output

Example:

  octave:1> for k = .1:.2:.5,
  > fprintf('1/%g = %10.2e\n',k,1/k); end
  1/0.1 =   1.00e+01
  1/0.3 =   3.33e+00
  1/0.5 =   2.00e+00

4 Graphics

4.1 2D graphics

Procedure for plotting a function $ y=f(x)$ :

  1. Generate a vector with the $ x$ -coordinates of the points to be plotted.
      x = x_min:step_size:x_max;
    
    (See also Section 2.1.)
  2. Generate a vector containing the corresponding $ y$ -values by letting $ f$ act on the vector $ x$ element-wise.
      y = f(x);
    
    Important: Since $ f$ operates element-wise, you must use of the operators .+, .-, .^ etc. instead of the usual +, - and ^! (See Section 2.4.)
  3. Finally you call the plot command.
      plot(x,y)
    
  4. To generate a coordinate grid, write:
      plot(x,y)
      grid
    

Example:

  octave:1> x = -10:.1:10;
  octave:2> y = sin(x).*exp(-abs(x));
  octave:3> plot(x,y)
  octave:4> grid
\includegraphics[height=7cm]{2d-plot1}

4.2 3D graphics:

Example:

  octave:1> x = -2:0.1:2;
  octave:2> [xx,yy] = meshgrid(x,x);
  octave:3> z = sin(xx.^2-yy.^2);
  octave:4> mesh(x,x,z);
\includegraphics[height=8cm]{3d-plot1}

4.3 Commands for 2D and 3D graphics

5 Exercises

5.1 Linear algebra

Take a matrix $ A$ and a vector $ b$ with

$ \displaystyle A = \begin{pmatrix}1 & 2  3 & 4
\end{pmatrix}$   and$\displaystyle \quad
b = \begin{pmatrix}
36 \\
88
\end{pmatrix}.
$

Solve the system of equations $ Ax = b$ . Calculate the LU and QR decompositions, and the eigenvalues and eigenvectors of $ A$ . Compute the Cholesky decomposition of $ A^T A$ , and verify that $ \operatorname{cond}(A^T A) = \operatorname{cond}(A)^2$ .

Solution:

  A = reshape(1:4,2,2).'; b = [36; 88];
  A\b
  [L,U,P] = lu(A)
  [Q,R] = qr(A)
  [V,D] = eig(A)
  A2 = A.'*A;
  R = chol(A2)
  cond(A)^2 - cond(A2)

5.2 Timing

Compute the matrix-vector product of a $ 100 \times 100$ random matrix with a random vector in two different ways. First, use the built-in matrix multiplication *. Next, use for-loops. Compare the results and computing times.

Solution:

  A = rand(100); b = rand(100,1);
  t = cputime;
  v = A*b; t1 = cputime-t;
  w = zeros(100,1);
  t = cputime;
  for n = 1:100,
    for m = 1:100
      w(n) = w(n)+A(n,m)*b(m);
    end
  end
  t2 = cputime-t;

  norm(v-w), t2/t1
Running this script yields the following output.
  ans = 0
  ans = 577.00

5.3 Stability functions of BDF-integrators

Calculate all the roots of the polynomial

$\displaystyle \frac{147}{60}\zeta^6 -6\zeta^5 + \frac{15}{2}\zeta^4
- \frac{20}{3}\zeta^3 + \frac{15}{4}\zeta^2 - \frac{6}{5}\zeta +
\frac{1}{6} .
$

Hint: Use the command compan.

Plot these roots as points in the complex plane and draw a unit circle for comparison. (Hint: hold, real, imag).

Solution:

  bdf6 = [147/60 -6 15/2 -20/3 15/4 -6/5 1/6];
  R = eig(compan(bdf6));
  plot(R,'+'); hold on
  plot(exp(pi*i*[0:.01:2]));
  if any(find(abs(R)>1))
    fprintf('BDF6 is unstable\n');
  else
    fprintf('BDF6 is stable\n');
  end

\includegraphics[height=8cm]{bdf6}

5.4 3D plot

Plot the graph of the function

$\displaystyle f(x,y) = \exp(-x^2-y^2) .
$

Solution:

  x = -3:0.1:3;
  [xx,yy] = meshgrid(x,x);
  z = exp(-xx.^2-yy.^2);
  figure, mesh(x,x,z);
  title('exp(-x^2-y^2)');

5.5 Hilbert matrix

For each $ n \times n$ Hilbert matrix $ H$ where $ n=1,\dots,15$ compute the solution to the linear system $ H x = b$ , $ b = $ ones(n,1). Calculate the error and the condition number of the matrix and plot both in semi-logarithmic coordinates. (Hint: hilb, invhilb.)

Solution:

  err = zeros(15,1); co = zeros(15,1);
  for k = 1:15
    H = hilb(k);
    b = ones(k,1);
    err(k) = norm(H\b-invhilb(k)*b);
    co(k) = cond(H);
  end
  semilogy(1:15,err,'r',1:15,co,'x');

5.6 Least square fit of a straight line

Calculate the least square fit of a straight line to the points $ (x_j,y_j)$ , given as two vectors $ x$ and $ y$ . Plot the points and the line.

Solution:

  function coeff = least_square (x,y)
    n = length(x);
    A = [x ones(n,1)];
    coeff = A\y;
    plot(x,y,'x');
    hold on
    interv = [min(x) max(x)];
    plot(interv,coeff(1)*interv+coeff(2));
  end

5.7 Trapezoidal rule

Write a program to integrate an arbitrary function $ f$ in one variable on an interval $ [a,b]$ numerically using the trapezoidal rule with $ h
= (b-a)/N$ :

$\displaystyle \int_a^b f(x)   dx \approx h\left( \frac{f(a)+f(b)}{2} +
\sum_{j=1}^{N-1} f(a+jh) \right).
$

For a function $ f$ of your choice, check, by generating a doubly logarithmic error plot, that the trapezoidal rule is of order $ 2$ .

Solution:

  function S = trapez(fun,a,b,N)
    h = (b-a)/N;
    % fy = feval(fun,[a:h:b]); better:
    fy = feval(fun,linspace(a,b,N+1));
    fy(1) = fy(1)/2;
    fy(N+1) = fy(N+1)/2;
    S = h*sum(fy);
  end

  function y = f(x)
    y = exp(x);
  end

  for k=1:15;
    err(k) = abs(exp(1)-1-trapez('f',0,1,2^k));
  end
  loglog(1./2.^[1:15],err);
  hold on;
  title('Trapezoidal rule, f(x) = exp(x)');
  xlabel('Increment');
  ylabel('Error');
  loglog(1./2.^[1:15],err,'x');