Octave Programming Tutorial/Print version

Octave Programming Tutorial

The current, editable version of this book is available in Wikibooks, the open-content textbooks collection, at
https://en.wikibooks.org/wiki/Octave_Programming_Tutorial

Permission is granted to copy, distribute, and/or modify this document under the terms of the Creative Commons Attribution-ShareAlike 3.0 License.

Getting started

The aim of this tutorial is to give you a quick introduction to basic Octave and to show that you know a lot of it already. If you should ever get stuck or need more information on an Octave function or command, type

help command

at the Octave prompt. command is the name of the Octave command or function on which to find help. Be warned though that the descriptions can sometimes be a bit technical and filled with jargon.

Starting Octave

Type octave in a terminal window to get started. You should see something like the following, depending on the version that you have installed:

GNU Octave, version 3.6.4.0
Copyright (C) 2013 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTIBILITY or
 FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/get-involved.html

Report bugs to <bug@octave.org> (but first, please read
http://www.octave.org/bugs.html to learn how to write a helpful report).

octave:1>

To exit Octave, type quit or exit.

Entering commands

The last line above is known as the Octave prompt and, much like the prompt in Linux, this is where you type Octave commands. To do simple arithmetic, use + (addition), - (subtraction), * (multiplication), / (division), and ^ (exponentiation).

Many mathematical functions are available and have obvious names (with many of those similar to other programming languages). For example, we have:

• trigonometric functions: sin, cos, tan
• inverse trigonometric functions: asin, acos, atan
• natural and base 10 logarithms: log, log10
• exponentiation: exp
• absolute value: abs

Various constants are also pre-defined: pi, e (Euler's number), i and j (the imaginary number sqrt(-1)), inf (Infinity), NaN (Not a Number - resulting from undefined operations, such as Inf/Inf.)

Here are some examples showing the input typed at the prompt and the output returned by Octave.

$2+3$	${\frac {\log _{10}100}{\log _{10}10}}$	$\left\lfloor {\frac {1+\tan 1.2}{1.2}}\right\rfloor$	${\sqrt {3^{2}+4^{2}}}$
octave:1> 2 + 3 ans = 5	octave:2> log10(100)/log10(10) ans = 2	octave:3> floor((1+tan(1.2)) / 1.2) ans = 2	octave:4> sqrt(3^2 + 4^2) ans = 5

$e^{i\cdot \pi }$

octave:5> e^(i*pi) ans = -1.0000e+00 + 1.2246e-16i octave:6> # Comment: From Euler's famous formula octave:6> # extremely close to the correct value of -1

Some things to note:

• Octave requires parentheses around the input of a function (so, log(10) is fine, but (log 10) is not).
• Any spacing before and after arithmetic operators is optional, but allowed.
• Not all Octave functions have obvious names (e.g. sqrt above). Don't panic for now. You will get to know them as we go along.

Plotting

You are going to plot the following pictures using Octave:

Figure 1	Figure 2
Figure 3	Figure 4

Figure 1 contains a plot of sin x vs x and is generated with the following commands. While this is a simple example, it is meant to illustrate the basic functionality. We will see more elaborate examples later.

x = linspace(0, 2*pi, 100);
y = sin(x);
plot(x, y);
figure;

The command that actually generates the plot is, of course, plot(x, y). Before executing this command, we need to set up the variables, x and y. The plot function simply takes two vectors of equal length as input, interprets the values in the first as x-coordinates and the second as y-coordinates and draws a line connecting these coordinates.

The first command above, x = linspace(0, 2*pi, 100), uses the linspace function to make a vector of equally spaced intervals (sometimes also called "linearly spaced values"). The first value in the vector is 0, the final value is 2π and the vector contains 100 values. This vector is assigned to the variable named x.

The second command computes the sin of each value in the vector variable, x, and stores the resulting vector in the variable y.

(As an aside: the name of a variable can be any sequence of letters, digits and underscores that does not start with a digit. There is no maximum length for variable names, and the case of alphabetical characters is important, i.e. a and A are two different variable names.)

Exercise

Plot the function $y=\lfloor x\rfloor$ for $x\in [0,15]$. (This is Figure 2).
Note : The graph may appear slightly inaccurate if the length(3rd) parameter of linspace is not sufficiently large. See the next heading for more information.

More on commands

The following commands and functions are useful for setting up variables for plotting 2D graphs.

• linspace creates a vector of evenly (linearly) spaced values.

Usage: linspace(start, stop, length). The length parameter is optional and specifies the number of values in the returned vector. If you leave out this parameter, the vector will contain 100 elements with start as the first value and stop as the last.

• plot plots a 2-dimensional graph.

Usage: plot(x, y) where x and y are vectors of equal length.

• figure creates a new plotting window.

This is useful for when you want to plot multiple graphs in separate windows rather than replacing your previous graph or plotting on the same axes.

• hold on and hold off sets whether you want successive plots to be drawn together on the same axes or to replace the previous plot.

• The ; at the end of a line suppress the display of the result. Try to remove it, and see !

Example

We are going to plot Figures 3 and 4. Figure 3 contains the 3 trigonometric functions

• $\cos 2x$,
• $\sin 4x$, and
• $2\sin x$

on one set of axes. Figure 4 contains the sum of these 3 functions.

Firstly, we use linspace to set up a vector of x-values.

octave:1> x = linspace(0, 2*pi);

Then, we compute the y-values of the 3 functions.

octave:2> a = cos(2*x);
octave:3> b = sin(4*x);
octave:4> c = 2*sin(x);

The following plots the first graph.

octave:5> figure;
octave:6> plot(x, a);
octave:7> hold on;
octave:8> plot(x, b);
octave:9> plot(x, c);

We use line 5 (figure) to tell Octave that we want to plot on a new set of axes. It is good practice to use figure before plotting any new graph. This prevents you from accidentally replacing a previous plot with the new one.

Note that on line 7, hold on is used to tell Octave that we don't want to replace the first plot (from line 6) with subsequent ones. Octave will plot everything after hold on on the same axes, until the hold off command is issued.

The figure you see shows all three plotted functions in the same color. To let Octave assign different colors automatically plot all functions in one step.

octave:10> plot(x, a, x, b, x, c);

Finally, we plot the second graph.

octave:11> figure;
octave:12> hold off;
octave:13> plot(x, a+b+c);

Line 11 creates a new graph window and line 12 tells Octave that any subsequent plots should simply replace previous ones. Line 13 generates the plot of the sum of the 3 trigonometric functions. A small note for mac users, if plot figure commands do not work first use this command: > setenv ("GNUTERM", "x11")

Exercises

• Plot the absolute value function for $x\in [-5,5]$.
• Plot a circle of radius 1, centered at the origin. (This is not so easy.)
• Plot your favourite function(s) like sin and cos

Warning

If you try (or have tried) to plot something like $x^{2}$ or $\sin x\times \cos x$, you will run into trouble. The following error messages are common. In the case of x^2:

error: for A^b, A must be square

In the case of sin(x)*cos(x):

error: operator *: nonconformant arguments (op1 is 1x100, op2 is 1x100)

This error occurs whenever you try multiply or divide two vector variables (remember that x and y are vectors). For now, you can do one of two things.

1. use plot(sin(x).*cos(x)) - notice the '.' between sin(x) and *. For x^2 plot(x.*x)
2. Read the section on vectors and matrices.

Challenge

Since Octave is a numerical (and not symbolic) mathematics package, it does make numerical errors and does not handle some operations well. To confirm this, make a plot of tan x, for x between -π and π. What is wrong with the resulting picture?

Your task is to generate the (much better looking) graph below using what you have learned so far and the axis function. axis can be used to adjust which part of the plot is actually displayed on screen. Use the command help axis to determine how this function works.

It might take some thinking to get the asymptote lines at $x=\pm \pi /2$ right. You can use help plot to find out how to plot dotted lines. (Try to make the plot yourself before looking at the solution below.)

Following commands could be used to generate the plot shown above.

octave:1> x_part_left  = linspace(-pi,  -pi/2-0.001, 100);
octave:2> x_part_mid   = linspace(-pi/2, pi/2, 100);
octave:3> x_part_right = linspace( pi/2+0.001, pi,   100);
octave:4> plot(x_part_left, tan(x_part_left));
octave:5> hold on;
octave:6> plot(x_part_mid, tan(x_part_mid));
octave:7> plot(x_part_right, tan(x_part_right));
octave:8> y_limit = 4;
octave:9> axis([-pi, pi, -y_limit, y_limit]);
octave:10> plot(linspace(-pi/2, -pi/2, 100), linspace(-y_limit, y_limit, 100), '.');
octave:11> plot(linspace( pi/2,  pi/2, 100), linspace(-y_limit, y_limit, 100), '.');
octave:12> hold off;

The horizontal plot range -π to π is split into three vectors such that singular points are skipped for the plot. In lines 4-7 the separate parts of the tan function are plotted. Thereafter, in line 8, we choose a limit for the y-axis and use it to constrain the vertical plot range (using the axis command [line 9]). Finally we add the asymptote lines at $x=\pm \pi /2$ in a dotted style (lines 10 & 11).

Script files

It is useful to be able to save Octave commands and return them later on. You might want to save your work or create code that can be reused (by yourself or somebody else). Such files are known as Octave script files. They should be saved with a .m extension so that Octave can recognise them. (The .m extension is used because MATLAB calls its script files M-files and Octave is based on MATLAB.)

To run an existing script in Octave, you have to be in the same directory as the script file and type in the name of the file without the .m in Octave. For example, if I have a script called myscript.m in an octave directory, the following two commands will execute the script.

chdir('~/octave'); % This changes to the octave directory
myscript;

Note that the chdir('~/octave') command is necessary only if you are not already inside that directory when running Octave.

In the following section you will be shown a number of new statements that you can use to make your Octave code much more powerful. A number of example script files are provided and you should save these into a directory for later use. A good idea is to create a directory called octave in your home directory and store all your Octave files in there.

---

Return to the Octave Programming tutorial index

Creating vectors and matrices

Here is how we specify a row vector in Octave:

octave:1> x = [1, 3, 2]
x =

  1  3  2

Note that

• the vector is enclosed in square brackets;
• each entry is separated by an optional comma. x = [1 3 2] results in the same row vector.

To specify a column vector, we simply replace the commas with semicolons:

octave:2> x = [1; 3; 2]
x =

  1
  3
  2

From this you can see that we use a comma to go to the next column of a vector (or matrix) and a semicolon to go to the next row. So, to specify a matrix, type in the rows (separating each entry with a comma) and use a semicolon to go to the next row.

octave:3> A = [1, 1, 2; 3, 5, 8; 13, 21, 34]
A =

   1   1   2
   3   5   8
  13  21  34

Operators

You can use the standard operators to

• add (+),
• subtract (-), and
• multiply (*)

matrices, vectors and scalars with one another. Note that the matrices need to have matching dimensions (inner dimensions in the case of multiplication) for these operators to work.

• The transpose operator is the single quote: '. To continue from the example in the previous section,

octave:4> A'
ans =

   1   3  13
   1   5  21
   2   8  34

(Note: this is actually the complex conjugate transpose operator, but for real matrices this is the same as the transpose. To compute the transpose of a complex matrix, use the dot transpose (.') operator.)

• The power operator (^) can also be used to compute real powers of square matrices.

Element operations

When you have two matrices of the same size, you can perform element by element operations on them. For example, the following divides each element of A by the corresponding element in B:

octave:1> A = [1, 6, 3; 2, 7, 4]
A =

  1  6  3
  2  7  4

octave:2> B = [2, 7, 2; 7, 3, 9]
B =

  2  7  2
  7  3  9

octave:3> A ./ B
ans =

  0.50000  0.85714  1.50000
  0.28571  2.33333  0.44444

Note that you use the dot divide (./) operator to perform element by element division. There are similar operators for multiplication (.*) and exponentiation (.^).

Let's introduce a scalar for future use.

a = 5;

The dot divide operators can also be used together with scalars in the following manner.

C = a ./ B

returns a matrix, C where each entry is defined by

$C_{ij}={\frac {a}{B_{ij}}}$

i.e. a is divided by each entry in B. Similarly

C = a .^ B

return a matrix with

$C_{ij}=a^{B_{ij}}$

Indexing

You can work with parts of matrices and vectors by indexing into them. You use a vector of integers to tell Octave which elements of a vector or matrix to use. For example, we create a vector

octave:1> x = [1.2, 5, 7.6, 3, 8]
x =

  1.2000  5.0000  7.6000  3.0000  8.0000

Now, to see the second element of x, type

octave:2> x(2)
ans = 5

You can also view a list of elements as follows.

octave:3> x([1, 3, 4])
ans =

  1.2000  7.6000  3.0000

This last command displays the 1st, 3rd and 4th elements of the vector x.

To select rows and columns from a matrix, we use the same principle. Let's define a matrix

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

  1  2  3
  4  5  6
  7  8  9

and select the 1st and 3rd rows and 2nd and 3rd columns:

octave:5> A([1, 3], [2, 3])
ans =

  2  3
  8  9

The colon operator (:) can be used to select all rows or columns from a matrix. So, to select all the elements from the 2nd row, type

octave:6> A(2, :)
ans =

  4  5  6

You can also use : like this to select all Matrix elements:

octave:7> A(:,:)
ans =
 
  1   2   3
  4   5   6
  7   8   9

Ranges

We can also select a range of rows or columns from a matrix. We specify a range with

start:step:stop

You can actually type ranges at the Octave prompt to see what the results are. For example,

octave:3> 1:3:10
ans =

   1   4   7  10

The first number displayed was start, the second was start + step, the third, start + (2 * step). And the last number was less than or equal to stop.

Often, you simply want the step size to be 1. In this case, you can leave out the step parameter and type

octave:4> 1:10
ans =

   1   2   3   4   5   6   7   8   9  10

As you can see, the result of a range command is simply a vector of integers. We can now use this to index into a vector or matrix. To select the $2\times 2$ submatrix at the top left of A, use

octave:4> A(1:2, 1:2)
ans =

  1  2
  4  5

Finally, there is a keyword called end that can be used when indexing into a matrix or vector. It refers to the last element in the row or column. For example, to see the last column in a Matrix, you can use

octave:5> A(:,end)
ans =

  3
  6
  9

Functions

The following functions can be used to create and manipulate matrices.

Creating matrices

• tril(A) returns the lower triangular part of A.

• triu(A) returns the upper triangular part of A.

• eye(n) returns the $n\times n$ identity matrix. You can also use eye(m, n) to return $m\times n$ rectangular identity matrices.

• ones(m, n) returns an $m\times n$ matrix filled with 1s. Similarly, ones(n) returns $n\times n$ square matrix.

• zeros(m, n) returns an $m\times n$ matrix filled with 0s. Similarly, zeros(n) returns $n\times n$ square matrix.

• rand(m, n) returns an $m\times n$ matrix filled with random elements drawn uniformly from $[0,1]$. Similarly, rand(n) returns $n\times n$ square matrix.

• randn(m, n) returns an $m\times n$ matrix filled with normally distributed random elements.

• randperm(n) returns a row vector containing a random permutation of the numbers $1,2,\ldots ,n$.

• diag(x) or diag(A). For a vector, x, this returns a square matrix with the elements of x on the diagonal and 0s everywhere else. For a matrix, A, this returns a vector containing the diagonal elements of A. For example,

octave:16> A = [1, 2, 3; 4, 5, 6; 7, 8, 9]
A =

  1  2  3
  4  5  6
  7  8  9

octave:17> x = diag(A)
ans =

  1
  5
  9

octave:18> diag(x)
ans =

  1  0  0
  0  5  0
  0  0  9

• linspace(a, b, n) returns a vector with n values, such that the first element equals a, the last element equals b and the difference between consecutive elements is constant. The last argument, n, is optional with default value 100.

octave:186> linspace(2, 4, 2)
ans =

  2   4

octave:187> linspace(2, 4, 4)
ans =

  2.0000   2.6667   3.3333   4.0000

octave:188> linspace(2, 4, 6)
ans =

  2.0000   2.4000   2.8000   3.2000   3.6000   4.0000

• logspace(a, b, n) returns a vector with n values, such that the first element equals $10^{a}$, the last element equals $10^{b}$ and the ratio between consecutive elements is constant. The last argument, n is optional with default value 50.

octave:189> logspace(2, 4, 2)
ans =

    100   10000

octave:190> logspace(2, 4, 4)
ans =

  1.0000e+02   4.6416e+02   2.1544e+03   1.0000e+04

octave:191> logspace(2, 4, 5)
ans =

  1.0000e+02   3.1623e+02   1.0000e+03   3.1623e+03   1.0000e+04

Other matrices

There are some more functions for creating special matrices. These are

• hankel (Hankel matrix),
• hilb (Hilbert matrix),
• invhilb (Inverse of a Hilbert matrix),
• sylvester_matrix (Sylvester matrix) - In v3.8.1 there is a warning: sylvester_matrix is obsolete and will be removed from a future version of Octave; please use hadamard(2^k) instead,
• toeplitz (Toeplitz matrix),
• vander (Vandermonde matrix).

Use help to find out more about how to use these functions.

Changing matrices

• fliplr(A) returns a copy of matrix A with the order of the columns reversed, for example,

octave:49> A = [1, 2, 3, 4; 5, 6, 7, 8; 9, 10, 11, 12] 
A =

   1   2   3   4
   5   6   7   8
   9  10  11  12

octave:50> fliplr(A)
ans =

   4   3   2   1
   8   7   6   5
  12  11  10   9

• flipud(A) returns a copy of matrix A with the order of the rows reversed, for example,

octave:51> flipud(A)
ans =

   9  10  11  12
   5   6   7   8
   1   2   3   4

• rot90(A, n) returns a copy of matrix A that has been rotated by (90n)° counterclockwise. The second argument, $n$, is optional with default value 1, and it may be negative.

octave:52> rot90(A)
ans =

   4   8  12
   3   7  11
   2   6  10
   1   5   9

• reshape(A, m, n) creates an $m\times n$ matrix with elements taken from A. The number of elements in A has to be equal to $mn$. The elements are taken from A in column major order, meaning that values in the first column ($A_{11},\ldots ,A_{m1}$) are read first, then the second column ($A_{12},\ldots ,A_{m2}$), etc.

octave:53> reshape(A, 2, 6)
ans =

   1   9   6   3  11   8
   5   2  10   7   4  12

• sort(x) returns a copy of the vector x with the elements sorted in increasing order.

octave:54> x = rand(1, 6)
x =

  0.25500  0.33525  0.26586  0.92658  0.68799  0.69682

octave:55> sort(x)
ans =

  0.25500  0.26586  0.33525  0.68799  0.69682  0.92658

Linear algebra

For a description of more operators and functions that can be used to manipulate vectors and matrices, find eigenvalues, etc., see the Linear algebra section.

---

Return to the Octave Programming Tutorial index

Plotting

Octave can work with gnuplot, Grace, PLplot. Some people deem PLplot is a replacement of the traditional gnuplot in Octave.

2D plots

plot(y)

If a single data argument is supplied, it is taken as the set of Y coordinates and the X coordinates are taken to be the indices of the elements, starting with 1.

plot(x, y)

• If the first argument is a vector and the second is a matrix, the vector is plotted versus the columns (or rows) of the matrix. (using whichever combination matches, with columns tried first.)
• If the first argument is a matrix and the second is a vector, the columns (or rows) of the matrix are plotted versus the vector. (using whichever combination matches, with columns tried first.)
• If both arguments are vectors, the elements of Y are plotted versus the elements of X.
• If both arguments are matrices, the columns of Y are plotted versus the columns of X. In this case, both matrices must have the same number of rows and columns and no attempt is made to transpose the arguments to make the number of rows match.

linspace(b, l, n)

results in a vector with n linearly spaced numbers between b and l. It is often used to produce the vector X in the two-argument form of plot.

Example

x = linspace(0, 10, 100);
y = sin(x);
plot(x, y);

Non-linear plots

semilogx, semiology, loglog, polar

Formatting

The plot may be formatted using an additional FMT argument:

plot (x, [y,] FMT, …)

In the argument any number of the following options may appear:

    `-'
          Set lines plot style (default).

    `.'
          Set dots plot style.

    `@'
          Set points plot style.

    `-@'
          Set linespoints plot style.

    `^'
          Set impulses plot style.

    `L'
          Set steps plot style.

    `N'
          Interpreted as the plot color if N is an integer in the range
          1 to 6.

    `NM'
          If NM is a two digit integer and M is an integer in the range
          1 to 6, M is interpreted as the point style.  This is only
          valid in combination with the `@' or `-@' specifiers.

    `C'
          If C is one of `"r"', `"g"', `"b"', `"m"', `"c"', or `"w"',
          it is interpreted as the plot color (red, green, blue,
          magenta, cyan, or white).

    `";title;"'
          Here `"title"' is the label for the key.

    `+'
    `*'
    `o'
    `x'
          Used in combination with the points or linespoints styles,
          set the point style.

     The color line styles have the following meanings on terminals that
     support color.

          Number  Gnuplot colors  (lines)points style
            1       red                   *
            2       green                 +
            3       blue                  o
            4       magenta               x
            5       cyan                house
            6       white           there exists

     The FMT argument can also be used to assign key titles.  To do so,
     include the desired title between semi-colons after the formatting
     sequence described above, e.g. "+3;Key Title;" Note that the last
     semi-colon is required and will generate an error if it is left
     out.

Furthermore, a set of property-value pairs may also be applied, simultaneously, to each "line" drawn by plot.

plot (x, [y,] FMT, property, value, …)

Example

Here instead of line plot we place a '+' marker on the graph for each datapoint, choose the first color (red), and add a legend.

x = linspace(0, 2*pi, 50);
y = sin(x);
plot(x, y, '+1;sin(x);');

The second example adds another line and some property-value pairs as well.

x = linspace(0, 2*pi, 100);
y = sin(x);
y2 = cos(x);
plot(x, y, '+1;sin(x);', "markersize", 10, x, y2, ";cos(x);", "markersize", 5, "marker", '*');

Commands related to plot

xlabel, ylabel, title and refresh

3D plots

mesh, meshgrid, surf

Contour plots

contour

contour (Z)

contour (Z, VN)

contour (X, Y, Z)

contour (X, Y, Z, VN)

contour (..., STYLE)

contour (H, ...)

[C, H] = contour (...)

    Plot level curves (contour lines) of the matrix Z, using the
    contour matrix C computed by `contourc' from the same arguments;
    see the latter for their interpretation.  The set of contour
    levels, C, is only returned if requested.  For example:
         x = 0:2;
         y = x;
         z = x' * y;
         contour (x, y, z, 2:3)
    The style to use for the plot can be defined with a line style
    STYLE in a similar manner to the line styles used with the `plot'
    command.  Any markers defined by STYLE are ignored.
    The optional input and output argument H allows an axis handle to
    be passed to `contour' and the handles to the contour objects to  be returned.

contourc

[C, LEV] = contourc (X, Y, Z, VN)

    Compute isolines (contour lines) of the matrix Z.  Parameters X, Y
    and VN are optional.
    The return value LEV is a vector of the contour levels.  The
    return value C is a 2 by N matrix containing the contour lines in
    the following format
         C = [lev1, x1, x2, ..., levn, x1, x2, ...
              len1, y1, y2, ..., lenn, y1, y2, ...]
    in which contour line N has a level (height) of LEVN and length of
    LENN.
    If X and Y are omitted they are taken as the row/column index of
    Z.  VN is either a scalar denoting the number of lines to compute
    or a vector containing the values of the lines.  If only one value
    is wanted, set `VN = [val, val]'; If VN is omitted it defaults to
    10.

    For example,
         x = 0:2;
         y = x;
         z = x' * y;
         contourc (x, y, z, 2:3)
              =>   2.0000   2.0000   1.0000   3.0000   1.5000   2.0000
              2.0000   1.0000   2.0000   2.0000   2.0000   1.5000

Images

colormap, image, imshow, hsv2rgb

>> help colormap 'colormap' is a function from the file C:\Octave\Octave-4.0.1\share\octave\4.0.1\ m\image\colormap.m

-- Function File: CMAP = colormap ()
-- Function File: CMAP = colormap (MAP)
-- Function File: CMAP = colormap ("default")
-- Function File: CMAP = colormap ("MAP_NAME")
-- Function File: CMAP = colormap (HAX, ...)
-- Command:  colormap MAP_NAME
-- Function File: CMAPS = colormap ("list")
-- Function File:  colormap ("register", "NAME")
-- Function File:  colormap ("unregister", "NAME")
    Query or set the current colormap.

    With no input arguments, `colormap' returns the current color map.

    `colormap (MAP)' sets the current colormap to MAP.  The colormap
    should be an N row by 3 column matrix.  The columns contain red,
    green, and blue intensities respectively.  All entries must be
    between 0 and 1 inclusive.  The new colormap is returned.

    `colormap ("default")' restores the default colormap (the `jet'
    map with 64 entries).  The default colormap is returned.

    The map may also be specified by a string, "MAP_NAME", where
    MAP_NAME is the name of a function that returns a colormap.

    If the first argument HAX is an axes handle, then the colormap for
    the parent figure of HAX is queried or set.

    For convenience, it is also possible to use this function with the
    command form, `colormap MAP_NAME'.

    `colormap ("list")' returns a cell array with all of the available
    colormaps.  The options "register" and "unregister" add or remove
    the colormap NAME from this list.

    See also: jet.

Saving and printing graphs

Print a graph to a printer or save it to a file

print("-Pprinter")
print("filename", "-ddevice")

Devices:

         `ps'
         `ps2'
         `psc'
         `psc2'
               Postscript (level 1 and 2, mono and color)

         `eps'
         `eps2'
         `epsc'
         `epsc2'
               Encapsulated postscript (level 1 and 2, mono and color)

         `ill'
         `aifm'
               Adobe Illustrator

         `cdr'
         `corel'
               CorelDraw

         `hpgl'
               HP plotter language

         `fig'
               XFig

         `dxf'
               AutoCAD

         `mf'
               Metafont

         `png'
               Portable network graphics

         `pbm'
               PBMplus
         `svg`
               Scalar Vector Graphics

If the device is omitted, it is inferred from the file extension, or if there is no filename it is sent to the printer as postscript.

See also

• Creating graphic file in Octave
• commons:Category:Images with Octave source code

---

Return to the Octave Programming tutorial index

The disp function

The disp function displays the value of a variable (scalar, vector, matrix, string, etc.) in the same way as simply typing the name of the variable does. For example,

octave:1> x = [1, 2, 3]
x =

  1  2  3

octave:2> disp(x)
  1  2  3

The name of the variable is, however, not displayed. You can also display the result of a computation without the ans = that normally precedes it.

octave:3> log(10)
ans = 2.3026
octave:4> disp(log(10))
2.3026

The output of disp depends on the format command.

octave:5> format long
octave:6> disp(log(10))
2.30258509299405

The displayed value can be printed to the screen, saved in a string or saved to a file using fdisp.

octave:7> s = disp(log(10))
s = 2.30258509299405

Note that s is a string containing the characters shown above.

File output

The fdisp function can be used to save values to a file. Before we can do this, we have to open a file. This is done using the fopen command.

• fopen(filename, mode) opens a file and returns an identifier for it. The filename argument is a string and can be the name of any new or existing file in the current directory. The mode argument is a string that specifies whether the file is opened for
• reading (r),
• writing (w), or
• appending (a).

When a file is opened for writing, its contents are erased and replaced with new data. To keep the existing data in a file and add to the end thereof, use the append mode.

octave:10> file_id = fopen('mydata.txt', 'w');

Here, file_id is simply the name of a variable that we use to tell Octave which file to write.

• fdisp(file_id, value) writes value to the file identified by file_id.

The output written to the file will appear exactly the same as the output from the disp command.

It is important to close a file after all data has been written to it. Closing the file tells Octave to finalise any output that might still be pending and frees up the file so that it can be opened by other users or programmes.

• fclose(file_id) closes the file identified by file_id.

The printf function

The printf function is considerably more powerful than disp and, consequently, a bit more complicated to use. With printf, you can define exactly what the output of a value should look like. This includes specifying

• the number of significant digits to display;
• the format of the number (integer, real, scientific, etc.);
• other output to display before or after the value.

Since there are so many different ways to format output using printf, we will discuss only the basics here using examples. For more complete information, type

doc printf

in Octave and page through the help using the spacebar key.

Outputting to screen, string or file

The printf function displays its output on the screen. Use the sprintf to return the result in a string and fprintf to write to a file. Note that the fprintf requires one additional parameter to specify the file identifier.

• printf(format, value, ...)
• sprintf(format, value, ...)
• fprintf(fid, format, value, ...)

Note that these functions can output more than one value at the same time--more on this in the next section.

The format string

Let's look at an example.

octave:18> x = 10.1;
octave:19> y = 5.5;
octave:20> z = 'test';
octave:21> printf('An integer: %i. A real: %f. This is a %s.\n', x, y, z);
An integer: 10. A real: 5.500000. This is a test.

The important part is the first argument to the printf function on line 21. It specifies what the output of printf should look like. Essentially, the percentage sign (%) indicates that a value should be placed at its position in the format string. In the format string

'An integer: %i. A real: %f. This is a %s.\n'

the %i is replaced with an integer, the %f with a real (f is for floating point) value, and the %s with a string. The values of the integer, real and string are given as arguments to printf after the format string. Note that x in the example above equals 10.1, but the value displayed is 10 since we specified that printf should display an integer. Finally, the \n at the end of the string tells Octave to move to a new line.

The next example demonstrates the following types:

• integer (%i),
• real (%f),
• scientific notation (%e),
• percentage symbol (%%).

For more types, see the Octave documentation.

octave:22> x = 10.34;
octave:23> printf("x is a real number: %f (%e in scientific notation).\n", x, x);
x is a real number: 10.340000 (1.034000e+01 in scientific notation).
octave:24> printf("Write percentages as %i%%.\n", x);
Write percentages as 10%.

Note that

• the %e format outputs a value in the form $a\times 10^{b}$, where $1\leq a<10$ and b is an integer;
• the variable x is passed to printf twice on line 23 since we want to output it twice (in different formats);
• the double percentage (%%) on line 24 outputs a single percentage symbol.

We can customize the output of values even further by specifying

• the width of the output, and
• the precision of the output.

The width allows you to right align numbers and is specified between the percentage and the format specifier. For example,

octave:36> x = 10;
octave:37> y = pi;
octave:38> z = 'test';
octave:39> printf("%9i\n%9f\n%9s\n", x, y, z);
       10
 3.141593
     test

Note that in the printf, each format specifier contains an integer (9). This tells printf to use 9 columns to output the integer, real and string.

The effect of the precision parameter depends on the type of output. It's most obvious use is to specify the number of digits to display after the decimal point of a real number. The precision is specified after the width and is preceded by a dot (.).

octave:40> printf("%9.3f\n", pi);
    3.142

This displays π using 9 columns and to 3 digits after the decimal point. Note that the number is rounded. For other uses of the precision parameter (e.g. for integers and strings), see the Octave help.

Outputting matrices

When the value passed to the printf function is a matrix or vector, all of the values in it are printed. If there is only one format specifier in the format string, it is used for each value in the matrix.

octave:51> A = [1, 2, 3; 4, 5, 6; 7, 8, 9];
octave:52> printf("%i\n", A)
1
4
7
2
5
8
3
6
9

Note that the values are read from the matrix in column-major order, i.e. all the values from the first column are displayed first, then the second column, etc..

If there is more than one format specifier in the format string, printf cycles through them.

octave:57> printf("[%i, %.1f, %.2e]\n", A)
[1, 4.0, 7.00e+00]
[2, 5.0, 8.00e+00]
[3, 6.0, 9.00e+00]

The values are still read in column-major order.

---

Return to the Octave Programming Tutorial index

General Mathematical Functions

Constants

• e is the base of the natural logarithm.

e without arguments returns the scalar e.

e(N) returns a square matrix of e of size N.

e(N, M, ...) where the arguments are dimensions of some matrix of e.

e(..., CLASS) where CLASS is an optional argument that specifies the return type, double or single.

• eps is the machine precision and returns the relative spacing between any floating point number and the next representable number. This value is system dependent.

eps returns the value of eps(1.0).

eps(X) returns the spacing between X and the next value.

eps with more than one argument is treated like e with the matrix value being eps(1.0).

• All of the constant functions listed are defined exactly like e

pi is the ratio of the circumference to the diameter of any circle.

I is the imaginary unit defined so I^2 = -1.

Inf is used for values that overflow the standard IEEE floating point range or the result of division by zero.

NaN is used for various results that are not well defined or undefined. Note that NaN never equals other NaN values. Use the function isnan to check for NaN.

realmax is the largest floating point value representable.

realmin is the smallest positive floating point value representable.

Arithmetic Functions

• floor(X) and ceil(X) return the highest integer not greater than X or lowest integer not less than X, respectively.
• round(X) and fix(X) return the integer closest to X or truncate X towards zero, respectively.
• rem(X,Y) and mod(X,Y) returns x - y * fix( x ./ y ) or x - y * floor( x ./ y ), they are the same except when dealing with negative arguments.
• hypot(X, Y) returns the length of the hypotenuse of a right-angle triangle with the adjacent and opposite of size X and Y.
• abs(X) return absolute of x.
• sign(X) return sign of the x (-1, 0 or +1).

Ordinary Trigonometry

• cos(X), sin(x) and tan(X) — the elemental functions that we all know and love. They take their arguments in radians.
• acos(X), asin(X) are the inverses of cos and sin and are able to compute arguments not contained in the range [-1,1].
• atan(X) and atan2(Y, X) are the 2 available inverses of tan. atan is a simple inverse whereas atan2 takes 2 arguments and returns an angle in the appropriate quadrant. More information on atan2 can be found here.
• Note that one can add the character d to any of the functions except atan2 and they will work in degrees rather than radians. For example: asind(0.3) = asin(0.3*180/pi)
• exp(x) , exponential function of x
• log(x) , natural logarithmic of x, log_e NOT log ₁₀

Hyperbolic Trigonometry

• cosh(X), sinh(X) and tanh(X) are analog to their more prosaic counterparts but deal with the unit hyperbola instead of the unit circle. They also take their arguments in radians.
• acosh(X), asinh(X) and atanh(X) are the inverses of cosh, sinh and tanh.
• Unlike their circular uncles they cannot be made to take their arguments in degrees.

Loops and conditions

Loops are used to repeat a block of code for a known or unknown number of times, depending on the type of loop. Using loops, you will draw some nice pictures of fractals and shapes drawn with random dots.

The for loop

We use for loops to repeat a block of code for a list of known values. As an example, we'll calculate the mean of a list of values. The mean is calculated from

${\overline {x}}={\frac {1}{n}}\sum _{i=1}^{n}x_{i}.$

We set up a vector with some values

octave:1> x = [1.2, 6.3, 7.8, 3.6];

and calculate the mean with

octave:2> sum = 0;
octave:3> for entry = x,
octave:4>   sum = sum + entry;
octave:5> end;
octave:6> x_mean = sum / length(x)

Line 2: Set sum equal to 0.
Line 3: For each value in x, assign it to entry.
Line 4: Increment sum by entry.
Line 5: Ends the for loop when there are no more members of x.
Line 6: Assign the final value of sum divided by the length of x to x_mean.

TO DO: get a better example and explain the code.

In general, we write a for loop as

for variable = vector
   ...
end

The ... represents the block of code that is executed exactly once for each value inside the vector.

Example: The Sierpinski triangle

The Sierpinski triangle is a fractal that can be generated with a very simple algorithm.

1. Start on a vertex of an equilateral triangle.
2. Select a vertex of the triangle at random.
3. Move to the point halfway between where you are now and the selected vertex.
4. Repeat from step 2.

Plotting the points that you visit by following this procedure, generates the following picture.

The code that generates this fractal is shown bellow. Note that this code uses one very simple for loop to generate the fractal (for i = 1:N ; ... ; end)

axis ([-1, 1, -0.75, 1.25], "square");
figure(1, "visible", "off");                          % no plotting window
hold on;

% defining the vertices of an equilateral triangle (symmetric to y-axis)
V = [ 0, 1;                                           % top vertex
     cos( (pi/2)+(2*pi/3) ), sin( (pi/2)+(2*pi/3) );  % left vertex
     cos( (pi/2)-(2*pi/3) ), sin( (pi/2)-(2*pi/3) )   % right vertex
     ];

r = floor(3*rand(1)+0.9999999999);                    % integer random number in [1:3]
x = [ V(r,1), V(r,2) ];                               % initializing x on random vertex

for i = 1:1000 % !!! 100000 => time=130m57.343s
  r = floor(3*rand(1)+0.9999999999);                  % integer random number in [1:3]
  x = ( x+V(r,[1:2]) )./2;                            % halfway, current to random vertex
  plot(x(1),x(2), ".");
endfor

print -dpng "sierpinski_m.png";

A Word of Caution

For performance reasons, it's best to perform as few tasks as possible inside 'for' loops, and graphical operations like plot should be moved outside of loops whenever possible. By simply storing all x values in a matrix and then plotting as shown below, the run time for the code to produce this figure drops to a few seconds, even when iterating over 100,000 elements (which the code above warns not to do).

axis ([-1, 1, -0.75, 1.25], "square");
figure(1, "visible", "off");                          % no plotting window
hold on;

% defining the vertices of an equilateral triangle (symmetric to y-axis)
V = [ 0, 1;                                           % top vertex
     cos( (pi/2)+(2*pi/3) ), sin( (pi/2)+(2*pi/3) );  % left vertex
     cos( (pi/2)-(2*pi/3) ), sin( (pi/2)-(2*pi/3) )   % right vertex
     ];

r = floor(3*rand(1)+0.9999999999);                    % integer random number in [1:3]
x(1,:) = [ V(r,1), V(r,2) ];                               % initializing x as a matrix this time.

for i = 2:100000 % Safe now: 100000 => time=7.85346s
  r = floor(3*rand(1)+0.9999999999);                  % integer random number in [1:3]
  x(i,:) = ( x(i-1,:)+V(r,[1:2]) )./2;                % Add each newly calculated x value to the matrix
endfor

plot(x(:,1),x(:,2), ".");                             % plot the entire matrix just once
print -dpng "sierpinski_m.png";

By initializing x as a matrix of the full final size ahead of time, this processing time can be reduced still further to only 1 or 2 seconds on modern hardware. In general, if a loop can be replaced with vectorization, it should be.

Exercises

1. Write a script that sums the first N integers. You can check your result with the formula ${\frac {1}{2}}N(N+1)$.
2. Write a script that does the same thing as the linspace function. It should start at some value, xstart, stop at xstop and create a vector that contains N values evenly spaced from xstart to xstop. You can use the zeros function to create a zero-filled vector of the right size. Use help zeros to find out how the function works.

The while loop

The while loop also executes a block of code more than once but stops based on a logical condition. For example

x = 1.0;
while x < 1000
   disp(x);
   x = x*2;
endwhile

will multiply x by 2 until its value exceeds 1000. Here, x < 1000 is the condition of the loop. As long as the condition holds (is true), the loop will continue executing. As soon as it is false, the loop terminates and the first instruction after the loop is executed.

The general form of a while loop is

while condition
   ...
endwhile

Exercise

1. Write a script that calculates the smallest positive integer, n, such that $a^{n}\geq b$ for some real numbers a and b. (Meaning, find the smallest power of a that is at least b.) Using the log function is considered cheating.

Example: The Mandelbrot fractal

The Mandelbrot set is another fractal and is generated by checking how long it takes a complex number to become large. For each complex number, c,

1. Start with $z_{0}=0$.
2. Let $z_{i}=z_{i-1}^{2}+c\quad \forall i=1,2,\ldots$
3. Find the first i such that $|z_{i}|>2$.

We record all of these i values and assign a colour to each of them. This is used to generate an image like this one.

You can download the code that generates this fractal from Mandelbrot.m. Note that there is a while loop (inside some for loops) that tests whether the complex number z has modulus less than 2:

while (count < maxcount) & (abs(z) < 2)
   ...
endwhile

The first condition in the while loop checks that we do not perform too many iterations. For some values of c the iteration will go on forever if we let it.

See also another version by Christopher Wellons

The do...until statement

These loops are very similar to while loops in that they keep executing based on whether a given condition is true or false. There are however some important difference between while and do...until loops.

1. while loops have their conditions at the beginning of the loop;
   do...until loops have theirs at the end.
2. while loops repeat as long as the condition is true;
   do...until loops continue as long as theirs is false.
3. while will execute 0 or more times (because the condition is at the beginning);
   do...until loops will execute 1 or more times (since the condition is at the end).

The general form of a do...until loop is

do
   ...
until condition

Exercise

Write a script that calculates the greatest common divisor (GCD) of two positive integers. You can do this using Euclid's algorithm.

Challenge

Write a script that generates random number pairs (a, b) that are distributed uniformly

1. over the disc $\left\{(x,y)|x^{2}+y^{2}\leq 1\right\}$ (the first image below);
2. as in the second image below

File:Octave uniform random circle.png File:Octave uniform random challenge.png

The break and continue statements

Sometimes it is necessary to stop a loop somewhere in the middle of its execution or to move on to the next value in a for loop without executing the rest of the loop code for the current value. This is where the break and continue statements are useful.

The following code demonstrates how the break statement works.

total = 0;
while true
   x = input('Value to add (enter 0 to stop): ');
   if x == 0
      break;
   endif
   total = total+x;
   disp(['Total: ', num2str(total)]);
endwhile

Without the break statement, the loop would keep executing forever since the condition of the while loop is always true. The break allows you to jump past the end of the loop (to the statement after the endwhile).

The break statement can be used in any loop: for, while or do...until.

The continue statement also jumps from the inside of a loop but returns to the beginning of the loop rather than going to the end. In a

1. for loop, the next value inside the vector will be assigned to the for variable (if there are any left) and the loop restarted with that value;
2. while loop, the condition at the beginning of the loop will be retested and the loop continued if it is still true;
3. do...until loop, the condition at the end of the loop will be tested and the loop continued from the beginning if it is still false.

As an example, the following code will fill the lower triangular part of a square matrix with 1s and the rest with 0s.

N = 5;
A = zeros(N); % Create an N x N matrix filled with 0s

for row = 1:N
   for column = 1:N
      if column > row
         continue;
      endif
      A(row, column) = 1;
   endfor
endfor

disp(A);

Note that the inner for skips (continues) over the code that assigns a 1 to an entry of A whenever the column index is greater than the row index.

The if statement

The general form of the if statement is

if condition1
   ...
elseif condition2
   ...
else
   ...
endif

If condition1 evaluates to true, the statements in the block immediately following the if are executed. If condition1 is false, the next condition (condition2 in the elseif) is checked and its statements executed if it is true. You can have as many elseif statements as you like. The final set of statements, after the else, is executed if all of the conditions evaluate to false. Note that the elseif and else parts of the if statement are optional.

The following are all valid if statements:

% Take the log of the absolute value of x
if x > 0
   y = log(x);
elseif x < 0
   y = log(-x);
else
   disp("Cannot take the log of zero.");
endif

x = input("Enter a value: ");
if x > 0
   disp("The number is positive");
endif

if x < 0
   disp("The number is negative");
endif

if x == 0
   disp("The number is zero");
endif

Example: The fractal fern

This algorithm is not quite complete. Have a look at the .m file available from http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4372&objectType=file.

The image to the right can be generated with the following algorithm:

  1. Let x1 and y1 be random values between 0 and 1.
  2. Choose one of the linear transformations below to calculate (xi+1, yi+1) from (xi, yi):
        1. xi+1 = 0
           yi+1 = 0.16yi
        2. xi+1 = 0.20xi − 0.26yi
           yi+1 = 0.23xi + 0.22yi + 1.6
        3. xi+1 = −0.15xi + 0.28yi
           yi+1 = 0.26xi + 0.24yi + 0.44
        4. xi+1 = 0.85xi + 0.04yi
           yi+1 = −0.04xi + 0.85yi + 1.6
     The first transformation is chosen if probability 0.01, the second and third with probability 0.07 each and the fourth with probability 0.85.
  3. Calculate these values for i up to at least 10,000.

You can download the code that generates this fractal as fracfern.m (this is disabled for now).

---

Return to the Octave Programming tutorial index

Writing functions

Syntax

In Octave, function definitions use the following syntax:

 function [return value 1, return value 2, ... ] = name( [arg1, arg2, ...] )
     body
 endfunction

Examples

Factorial Function

The factorial function, which takes exactly one argument and returns one integer, is as follows.

 function result = factorial( n )
     if( n == 0 )
         result = 1;
         return;
     else
         result = prod( 1:n );
     endif
 endfunction

Save m-File

You can save the definition of the function in your working directory. Use the name of the function also for the filename - e.g. use the filename factorial.m for the definition of the function factorial .

Maximum and Minimum of two integer Values

The following function, maxmin, returns the maximum and minimum value of two integers:

 function [max,min] = maxmin( a, b )
    if(a >= b )
         max = a;
         min = b;
         return;
    else
         max = b;
         min = a;
         return;
    endif
 endfunction

To call a function with multiple arguments, you specify multiple variables to hold the results. For example, one could write:

 [big,small] = maxmin(7,10);

After executing, the variable 'big' would store the value '10' and the variable 'small' would store '7'. If fewer than two variables are used to hold the result then fewer outputs are returned. Writing

a = maxmin(13,5)

would store 13 in the variable 'a', and would discard the value of 5.

---

Return to the Octave Programming tutorial index

Vectorization

Writing routines in terms of vector operations can be orders of magnitudes more efficient than using built-in interpreted loops because Octave can make use of highly optimized FORTRAN and C numerical linear algebra libraries instead. Even if a routine or function is not written in a vectorized form, it is possible to take advantage of vectorization by using arrayfun or a similar structure.

vectorizing a regular function with arrayfun

Consider an anonymous function

octave:1> f = @(x) sin(x)*x

Octave output :

f =
@(x) sin (x)*x

and assume that we want to calculate this function for every element of a given vector of integers from 1 to 7 :

octave:2> y=1:7

y =
  1   2   3   4   5   6   7

then passing y as an argument for f will give error

octave:3> f(y)
error: operator *: nonconformant arguments (op1 is 1x7, op2 is 1x7)
error: called from:
error:    at line -1, column -1

this is because f is not defined for a vector input. But this is not a problem as we can do:

octave:4> arrayfun(f,y)

and output is :

ans =
  0.84147   1.81859   0.42336  -3.02721  -4.79462  -1.67649   4.59891

This is an order of magnitude faster than calculating f for many y values by using a loop which has a big overhead.

Functions

• d = det(A) computes the determinant of the matrix A.
• lambda = eig(A) returns the eigenvalues of A in the vector lambda, and
• [V, lambda] = eig(A) also returns the eigenvectors in V but lambda is now a matrix whose diagonals contain the eigenvalues. This relationship holds true (within round off errors) A = V*lambda*inv(V).
• inv(A) computes the inverse of non-singular matrix A. Note that calculating the inverse is often 'not' necessary. See the next two operators as examples. Note that in theory A*inv(A) should return the identity matrix, but in practice, there may be some round off errors so the result may not be exact.
• A / B computes X such that $XB=A$. This is called right division and is done without forming the inverse of B.
• A \ B computes X such that $AX=B$. This is called left division and is done without forming the inverse of A.
• norm(A, p) computes the p-norm of the matrix (or vector) A. The second argument is optional with default value $p=2$.
• rank(A) computes the (numerical) rank of a matrix.
• trace(A) computes the trace (sum of the diagonal elements) of A.
• expm(A) computes the matrix exponential of a square matrix. This is defined as

$I+A+{\frac {A^{2}}{2!}}+{\frac {A^{3}}{3!}}+\cdots$

• logm(A) computes the matrix logarithm of a square matrix.
• sqrtm(A) computes the matrix square root of a square matrix.

---

Below are some more linear algebra functions. Use help to find out more about them.

• balance (eigenvalue balancing),
• cond (condition number),
• dmult (computes diag(x) * A efficiently),
• dot (dot product),
• givens (Givens rotation),
• kron (Kronecker product),
• null (orthonormal basis of the null space),
• orth (orthonormal basis of the range space),
• pinv (pseudoinverse),
• syl (solves the Sylvester equation).

Factorizations

• R = chol(A) computes the Cholesky factorization of the symmetric positive definite matrix A, i.e. the upper triangular matrix R such that $R^{T}R=A$.
• [L, U] = lu(A) computes the LU decomposition of A, i.e. L is lower triangular, U upper triangular and $A=LU$.
• [Q, R] = qr(A) computes the QR decomposition of A, i.e. Q is orthogonal, R is upper triangular and $A=QR$.

---

Below are some more available factorizations. Use help to find out more about them.

• qz (generalized eigenvalue problem: QZ decomposition),
• qzhess (Hessenberg-triangular decomposition),
• schur (Schur decomposition),
• svd (singular value decomposition),
• housh (Householder reflections),
• krylov (Orthogonal basis of block Krylov subspace).

---

Return to the Octave Programming Tutorial index

Differential equations

With Octave there are two functions available, that can be used for solving differential equations. The implementation is intergated in the core functionality of Octave. The used algorithms are based on Solvers for Ordinary Differential Equations written in Fortran.

Ordinary Differential Equations

The following function lsode can be used for Ordinary Differential Equations (ODE) of the form using Hindmarsh's ODE solver LSODE.

Function: lsode (fcn, x0, t_out, t_crit)

• The first argument is the name of the function to call to compute the vector of right hand sides. It must have the form

   x_dot = f (x, t)

in this context x_dot is the return value of the function $f$ and $x$ is a vector (in the example below it has 2 components) and $t$ is a scalar. The solver lsode gets a vector of scalars as input.

• The second argument x0 defines the initial condition.
• By the third argument the output times t_out are defined by a vector. At these points in time the solution is requested for calculation. Points in time include as first element the corresponding time for the initial state.

• The last argument is optional parameter, that may be used to determine a set of points in times that the ODE solver should not integrate for. Especially when the solution contain singularities, these parameter might be used for successful run of the solver. Other use-case might be a discontinuity of the derivative.

Example

The following example can be used for solving a set of two differential equations using lsode. The function is defined by:

   function x_dot = f (x, t) 
 
     r = 0.15;
     k = 1.6;
     a = 1.25;
     b = 0.12;
     c = 0.89;
     d = 0.58;
 
     x_dot(1) = r*x(1)*(1 - x(1)/k) - a*x(1)*x(2)/(1 + b*x(1));
     x_dot(2) = c*a*x(1)*x(2)/(1 + b*x(1)) - d*x(2);
  
   endfunction

the function f is integrated with the command

   x = lsode ("f", [1; 2], (t = linspace (0, 50, 200)'));

The linspace command is producing a set of 200 scalar values stored in the variable t.

   plot (t, x)

Partial Differential Equations

For Partial Differential Equations (PDE) see the following PDE example.

John Weatherwax (2006) provided the octave code defined the derivatives, boundary conditions with m-files:

• DBDU_BNDRY.m,
• DBDUX_BNDRY.m,
• DZDT_BNDRY.m and
• UINIT.m with the initial system for the PDE.

Display the libraries with the links above and adapt the code to your PDE problem.

• pdecol_Script.m call the m-files with the boundary conditions. Defining the files in that way encapsules the definition of functions in separate m-files for the solver script pdecol_Script.m. Apply the script to a new problem requires the adaptation of problem specific m-files above and the solver script remains untouched.

Polynomials

In Octave, a polynomial is represented by its coefficients (arranged in descending order). For example, the vector

octave:1> p = [-2, -1, 0, 1, 2];

represents the polynomial

$-2x^{4}-x^{3}+x+2.$

You can check this by displaying the polynomial with the function polyout.

octave:2> polyout(p, 'x')
-2*x^4 - 1*x^3 + 0*x^2 + 1*x^1 + 2

The function displays the polynomial in the variable specified (x in this case). Note that the ^ means raised to the power of much like the Octave operator.

Functions

• Evaluating a polynomial:

y = polyval(p, x)

This returns $p(x)$. If x is a vector or matrix, the polynomial is evaluated at each of the elements of x.

• Multiplication:

r = conv(p, q)

Here, p and q are vectors containing the coefficients of two polynomials and the result, r, will contain the coefficients of their product.

(As an aside, the reason this function is called conv is that we use vector convolution to do polynomial multiplication.)

• Division:

[b, r] = deconv(y, a)

This returns the coefficients of the polynomials b and r such that

$y=ab+r.$

So, b contains the coefficients of the quotient and r the coefficients of the remainder of y and a .

• Root finding:

roots(p)

This returns a vector containing all the roots of the polynomial with coefficients in p.

• Derivative:

q = polyder(p)

This returns the coefficients of the derivative of the polynomial whose coefficients are given by vector p.

• Integration:

q = polyint(p)

This returns the coefficients of the integral of the polynomial whose coefficients are represented by the vector p. The constant of integration is set to 0.

• Data fitting:

p = polyfit(x, y, n)

This returns the coefficients of a polynomial $p(x)$ of degree $n$ that best fits the data $(x_{i},y_{i})$ in the least squares sense.

Warning: Adding polynomials

Since polynomials are represented by vectors of their coefficients, adding polynomials is not straightforward in Octave. For example, we define the polynomials $p(x)=x^{2}-1$ and $q(x)=x+1$.

octave:1> p = [1, 0, -1];
octave:2> q = [1, 1];

If we try to add them, we get

octave:3> p + q
error: operator +: nonconformant arguments (op1 is 1x3, op2 is 1x2)
error: evaluating binary operator `+' near line 22, column 3

This happens because Octave is trying to add two vectors (p and q) of different lengths. This operation is not defined. To work around this, you have to add some leading zeroes to q.

octave:4> q = [0, 1, 1];

Note that adding leading zeroes does not change the polynomial. Next, we add

octave:5> p + q
ans =

  1  1  0

octave:6> polyout(ans, 'x')
1*x^2 + 1*x^1 + 0

---

Return to the Octave Programming tutorial index

Sets

SET

Octave set operation There are different set operations in octave basically Octave can use vector, cell arrays or matrix as its set input operation.

1) SET UNION:

The set union operations is one of the operations where two sets a and b are merged together.

In this example we will took two set namely ,A with contents 1,2,3 and B with contents 3,4,5. First we are going to use the set union operation using Octave. This can be accomplished using the union function available with octave.

 octave:1> a=[1,2,3]  
 a =  
   1  2  3  
 octave:2> b=[3,4,5]  
 b =  
   3  4  5  
 octave:3> union(a,b)  
 ans =  
   1  2  3  4  5

2) SET INTERSECTION :

The set intersection operations is one of the operations where two sets a and b are merged together and common element is taken.

In this example we will took two set namely A with contents 1,2,3 and b with contents 3,4,5 .First we are going to use the set intersection operation using Octave . This can be accomplished using the intersection function available with octave.

octave:1> a=[1,2,3]  
 a =  
   1  2  3  
 octave:2> b=[3,4,5]  
 b =  
   3  4  5  
 octave:3> intersect(a,b)  
 ans =  
   3

3) SET DIFFERENCE:

The set difference operations also called as the a-b operation is the operation which returns those element of a that are not in b. Lets write down the set difference operation as fallows :

Set=A-B

In this example we will took two set namely A with contents 1,2,3 and b with contents 3,4,5 First we are going to use the setdifference operation using Octave . This can be accomplished using the setdiff function available with octave. The difference operation is a Set operation in which those elements of a that are not in b are returned .

 octave:1> a=[1,2,3]  
 a =  
   1  2  3  
 octave:2> b=[3,4,5]  
 b =  
   3  4  5  
 octave:3> setdiff(a,b)  
 ans =  
   1,2

4) UNIQUE:

The unique operation returns a set containing one copy of each element in the given set, and with the elements in sorted order.

In this example we will take a set namely A with contents 1,2,4,3,5,2,1 containing the duplicate elements 1 and 2, and elements out of order. The unique operation will sort the set and remove extra copies of the duplicated elements.

 octave:1> a=[1,2,4,3,5,2,1]  
 a =
   1   2   4   3   5   2   1  
 octave:2> unique(a)
 ans =
   1   2   3   4   5
</syntaxhihglight>


4) XOR operation:
The Xor operation is one of the operation in which the two sets are compared  and if both the set evalulation remains the same then the result is return false else its return true . To execute set exor operation in octave we are going to use setxor function.
<syntaxhighlight lang="octave">
  a=[1,2,3]  
  a =  
    1  2  3  
  octave-3.2.4:34> b=[1,2,4]  
  b =  
    1  2  4  
  octave-3.2.4:35> setxor(a,b)  
  ans =  
    3  4

5) Ismember

The octave set ISMEMBER is one the function in which two sets are compared and the those elements that are present in the second set are marked as true rest are marked as false. This function is basically used to check which element are present in both the set.

 ocatave-3.2.4:33> a=[1,2,3]  
 a =  
   1  2  3  
  octave-3.2.4:34> b=[1,2,4]  
  b =  
    1  2  4  
  octave-3.2.4:35>  ismember(a,b)  
  ans =  
    [1 1 0]