# Octave Programming Tutorial/Print version

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

# 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
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:
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.
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.
Then, we compute the y-values of the 3 functions.
The following plots the first graph.
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.
Finally, we plot the second graph.
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:
In the case of sin(x)*cos(x):
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.
* use plot(sin(x).*cos(x)) - notice the '.' between sin(x) and *. For x^2 plot(x.*x)
* 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.
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.
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

## References

* ↑ As usual in Mathematics, the trigonometric functions take their arguments in radians.

# Creating vectors and matrices

Here is how we specify a row vector in Octave:
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:
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.

# 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,
(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:
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.
The dot divide operators can also be used together with scalars in the following manner.
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
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
Now, to see the second element of x, type
You can also view a list of elements as follows.
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
and select the 1st and 3rd rows and 2nd and 3rd columns:
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
You can also use : like this to select all Matrix elements:

## Ranges

We can also select a range of rows or columns from a matrix. We specify a range with
You can actually type ranges at the Octave prompt to see what the results are. For example,
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
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
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

# 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,
* 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.
* 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.

### 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,
* flipud(A) returns a copy of matrix A with the order of the rows reversed, for example,
* 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.
* 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.
* sort(x) returns a copy of the vector x with the elements sorted in increasing order.

# 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

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.
* 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.
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.

### Non-linear plots

semilogx, semiology, loglog, polar

### Formatting

The plot may be formatted using an additional FMT argument:
In the argument any number of the following options may appear:
Furthermore, a set of property-value pairs may also be applied, simultaneously, to each "line" drawn by plot.

### 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.
The second example adds another line and some property-value pairs as well.

### 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 (...)
contourc
[C, LEV] = contourc (X, Y, Z, VN)

## 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

## Saving and printing graphs

Print a graph to a printer or save it to a file
Devices:
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,
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.
The output of disp depends on the format command.
The displayed value can be printed to the screen, saved in a string or saved to a file using fdisp.
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.
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
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.
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
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.
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,
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 (.).
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.
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.
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 $_{10}$

## 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
and calculate the mean with
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
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.
* Start on a vertex of an equilateral triangle.
* Select a vertex of the triangle at random.
* Move to the point halfway between where you are now and the selected vertex.
* 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)

#### 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).
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

* Write a script that sums the first N integers. You can check your result with the formula ${\frac {1}{2}}N(N+1)$.
* 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
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

### Exercise

* 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,
* Start with $z_{0}=0$.
* Let $z_{i}=z_{i-1}^{2}+c\quad \forall i=1,2,\ldots$
* 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:
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.
* while loops have their conditions at the beginning of the loop; do...until loops have theirs at the end.
* while loops repeat as long as the condition is true; do...until loops continue as long as theirs is false.
* 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

### 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
* over the disc $\left\{(x,y)|x^{2}+y^{2}\leq 1\right\}$ (the first image below);
* as in the second image below

## 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.
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
* 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;
* while loop, the condition at the beginning of the loop will be retested and the loop continued if it is still true;
* 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.
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 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:

### 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:
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:

## Examples
### Factorial Function

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

### 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:
To call a function with multiple arguments, you specify multiple variables to hold the results. For example, one could write:
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
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 output :
and assume that we want to calculate this function for every element of a given vector of integers from 1 to 7 :
then passing y as an argument for f will give error
this is because f is not defined for a vector input. But this is not a problem as we can do:
and output is :
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
    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:
the function f is integrated with the command
The linspace command is producing a set of 200 scalar values stored in the variable t.

### Options for lsode
## 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.

## References
# Polynomials

In Octave, a polynomial is represented by its coefficients (arranged in descending order). For example, the vector
represents the polynomial
$-2x^{4}-x^{3}+x+2.$
You can check this by displaying the polynomial with the function polyout.
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:
This returns $p(x)$. If x is a vector or matrix, the polynomial is evaluated at each of the elements of x.
* Multiplication:
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:
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:
This returns a vector containing all the roots of the polynomial with coefficients in p.
* Derivative:
This returns the coefficients of the derivative of the polynomial whose coefficients are given by vector p.
* Integration:
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:
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$.
If we try to add them, we get
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.
Note that adding leading zeroes does not change the polynomial. Next, we add
Return to the Octave Programming tutorial index

# Sets

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.
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.
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 .
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.
* ↑ Alan C. Hindmarsh (1983), ODEPACK, A Systematized Collection of ODE Solvers, in Scientific Computing, R. S. Stepleman, editor,
* ↑ John Weatherwax (2006) Solving Partial Differential Equations in Octave URL: https://waxworksmath.com/Software/PDECOL/Web/pdecol_example1.html (accessed 2020/07/14)