3
Scientific Programming

This text contains basic physics programs written in the Octave scientific programming language that is freely available from http://www.gnu.org/software/octave/index.html with documentation at www.octave.org. Default selections can be chosen during setup. Octave incorporates many features of the commercial MATLAB® language and facilitates rapid and compact coding (for a more extensive introduction, refer to A Short Course in Computational Science and Engineering: C++, Java and Octave Numerical Programming with Free Software Tools, by David Yevick Copyright © 2012 David Yevick). Some of the material in the following text is reprinted with permission from Cambridge University Press.

3.1 Language Fundamentals

A few important general programming concepts as applied to Octave are first summarized below:

1. A program consists primarily of statements that result from terminating a valid expression not followed by the continuation character . (three lower dots), a carriage return, or a semicolon.
2. An expression can be formed from one or more subexpressions linked by operators such as + or *.
3. Operators possess different levels of precedence, e.g., in 2/4 + 3, the division operation possesses a higher precedence and is therefore evaluated before addition. In expressions involving two or more operators with the same precedence level, such as division and multiplication, the operations are typically evaluated from left to right, e.g., 2/4 * 3 equals (2/4) * 3.
4. The parenthesis operator, which evaluates the expression that it encloses, is assigned to the highest precedence level. This eliminates errors generated by incorrect use of precedence or associativity.
5. Certain style conventions, while not required, enhance clarity and readability:
   1. Variables and function names should be composed of one or more descriptive words. The initial letter should be uncapitalized, while the first letter of each subsequent word should be capitalized as in outputVelocity.
   2. Spaces should be placed to the right and left of binary operators, which act on the expressions (operands) to their left and right, as in 3 + 4, but no space should be employed in unary operator such as the negative sign in -3 + 4. Spaces are preferentially be inserted after commas as in computeVelocity( 3, 4 ) and within parentheses except where these indicate indices.
   3. Indentation should be employed to indicate when a group of inner statements is under the logical control of an outer statement such as in if ( firstVariable == 0 ) secondVariable = 5; end
   4. Any part of a line located to the right of the symbol % constitutes a comment that typically documents the program. Statements that form a logical unit should be preceded by one or more comment lines and surrounded by blank lines. Statement lines that introduce input variables should end with a comment describing the variables.

3.1.1 Octave Programming

Running Octave: Starting Octave opens a command window into which statements can be entered interactively. Alternatively, a program in the directory programs in partition C: is created by first entering cd C:\programs into the command window, pressing the enter key, and then entering the command edit. Statements are then typed into the program editor, the file is saved by selecting Save from the button or menu bar as a MATrix LABoratory file such as myFile.m (the .m extension is appended automatically by the editor), and the program is then run by typing myFile into the command window. The program can also be activated by including the statement myFile; within another program. To list the files in the current directory, enter dir into the Octave command window.

Help Commands: Typing help commandName yields a description of the command commandName. To find all commands related to a word subject, typelookfor subject. Entering doc or doc topic brings up, respectively, a complete help document and a description of the language feature topic.

Input and Output: A value of a variable G can be entered into a program (.m file) from the keyboard by including the line G = input( 'user prompt' ). The statement format long e sets the output style to display all 15 floating-point number significant digits, after which format short e reverts to the default 5 output digits.

Constants and Complex Numbers: Some important constants are i and j, which both equal , e, and pi. However, if a variable assignment such as i = 3; is encountered in an Octave program, i ceases to be identified with the imaginary unit until the command clear i is issued. Imaginary numbers can be manipulated with the functions real( ), imag( ), conj( ), and norm( ), and imaginary values are automatically returned by standard functions such as exp( ), sin( ), and sinh( ) for imaginary arguments.

Arrays and Matrices: A symbol A can represent a scalar, row, or column vector or matrix of any dimension. Row vectors are constructed either by

or

The corresponding column vector can similarly be entered in any of the following three ways:

Here .' indicates transpose, while ' instead implements the Hermitian (complex conjugate) transpose.

A 2 × 2 matrix

can be constructed by, e.g., mRC = [ 1 2; 3 4 ]; after which mRC(1, 2) returns (MRC)12, here the value 2. Subsequently, size(mRC) yields a vector containing the row and column dimensions of mRC, while length( mRC ) returns the maximum of these values. Here, we introduce the convention of appending R, C, or RC to the variable name to respectively identify row vectors, column vectors, and matrices.

Basic Manipulations: A value n is raised to the power m by n^m. The remainder of n/m is denoted rem( n, m ) and is positive or zero for n > 0 and negative or zero for n < 0. The function mod( n, m ) returns n modulus m, which is always positive, while ceil( ), floor( ), and fix( ) round floating-point numbers to the next larger integer, smaller integer, and nearest integer closer to zero, respectively.

Vector and Matrix Operations: Two vectors or matrices of the same dimension can be added or subtracted. Multiplying a matrix or vector by a scalar, c, multiplies each element by c. Additionally, eye( n, n ) is the n × n unit or identity matrix with ones along the main diagonal and zeros elsewhere, while ones( n, m ) and zeros( n, m ) are n × m matrices with all elements one or zeros so that

and

Further, mRC * mRC, or equivalently mRC^2, multiplies mRC by itself, while

implements component-by-component multiplication. Other arithmetic operations function analogously so that the (i, j) element of M ./ N is M ij /N ij . Functions such as cos( M ) return a matrix composed of the cosines of each element in M.

Solving Linear Equation Systems: The solution of the linear equation system xR * mRC = yR is xR = yR / mRC, while mRC * xC = yC is solved by xC = mRC \ yC. The inverse of a matrix mRC is represented by inv( mRC ). The eigenvalues of a matrix are obtained through eigenValues = eig( mRC ), while both the eigenvalues and eigenvectors are returned through [ eigenValues, eigenVectors ] = eig( mRC ).

Random Number Generation: A single random number between 0 and 1 is generated by rand, while rand( m, n ) returns a m × n matrix with random entries. The same random sequence can be generated each time a program is run by including rand( 'state', 0 ) before the first call to rand.

Control Logic and Iteration: The logical operators in octave are ==, <, <=, >, >=, ~= (not equal) and the and, or, and not operators-&, |, and ~, respectively. Any nonzero value is taken to represent a logical "true" value, while a zero value corresponds to a logical "false" as can be seen by evaluating, e.g., 3 & 4, which produces the output 1. Thus,

executes the statements denoted by xxx if the logical statement S == 2 is true, yyy if S == 3, and zzz otherwise. The for loop

yields the array vR = [ sin( p ) sin( 9p / 10 ) . sin( p/10 ) 0 ], while 1 : 10 yields an array with elements from 1 to 10 in unit increments. Mistakenly replacing colons by commas or semicolons results in severe and often difficult to detect errors. If a break statement is encountered within...