Wednesday, May 7, 2008

Print to png file from Octave

His last advice is very useful. I tried and it works perfectly.
Here is the link to the original post.

----------------------------------------------------------------

I ran into a problem trying to print a plot to a PNG file. I originally reported this to the octave bugs list but I finally realized it was more a matter of understanding how to use fonts with Octave 3.0 under Linux. (I don't know if this problem exists for Windows users or if so how to solve it). It took me a long time to figure out how this works so I thought I'd share the answer.

I was trying to do something like this:

x = 1:10 ;
plot (x);
title ('My plot');
print myplot.png -dpng

When Octave executed the print command, I got this error message:

gdImageStringFT: Could not find/open font while printing string My plot with font Helvetica

It turns out that Octave 3.0 uses 'Helvetica' as it's default font and if you don't happen to have a Helvetica font installed, you'll get this error message. I fixed this problem by finding the fonts that are installed on my system and explicitly specifying which font to use. Here's a detailed explanation:

First you need some fonts. You could use one of the five default fonts available in gnuplot but I found them way too limited for my needs. You probably already have some good true type fonts installed in your Linux distro. If not, you should be able to install some with your package manager.

Next you need to find the path to the fonts. Try looking in "/usr/share/fonts/truetype". If they aren't there, you can do "find . -name *.ttf" from /usr or other system directories until you find the font directories. As an example, I found two font files I wanted to use:

/usr/share/fonts/truetype/ttf-mgopen/MgOpenModernaBold.ttf
/usr/share/fonts/truetype/ttf-dejavu/DejaVuSans-Bold.ttf

Now you need to tell gnuplot where to find the fonts. Do this by setting the GDFONTPATH environment variable. I had to list each font subdirectory like this:

export GDFONTPATH="/usr/share/fonts/truetype/ttf-mgopen:/usr/share/fonts/truetype/ttf-dejavu"

Now you can use the name of the font file in Octave like this:

title ('My plot', 'FontName', 'MgOpenModernaBold.ttf', 'FontSize', 18);

As far as I can tell, the Octave plot window didn't use the font I specified but when I printed to a PNG file, it had the right font.

If for some reason, setting GDFONTPATH doesn't work, you can specify the full path to the font. It's a yucky hack, but it works. I discovered thru trial and error that there seems to be some magical length limit on the path name. When I did this in octave:

title ('My plot', 'FontName', "/usr/share/fonts/truetype/ttf-mgopen/MgOpenModernaBold.ttf");

I got an error. So I copied the font file to a directory with a shorter name and shortened the file name as well:

title ('My plot', 'FontName', "/home/dave/Fonts/MdrnaBd.ttf");

This works just fine.

Monday, May 5, 2008

Octave Tutorial

Link

http://homepages.nyu.edu/~kpl2/dsts6/octaveTutorial.html

Octave resources:

Type help help at the octave command line to find out about the built in documentation system

Type help > to get documentation on a particular keyword, operator, or function. Here is an example of asking for documentation on the sin() function: octave:1>help sin

On-line Tutorial: There are several on-line octave/MATLAB tutorials. Here is a link to one which gives a good introduction to the language: http://www.aims.ac.za/resources/tutorials/octave/

On-line octave Manual: The octave manual is available as a PDF here: http://www.octave.org/docs.html

octave-forge documentation: octave-forge is a library of octave functions for working on domain specific (including audio signal processing) problems. You can install octave-forge using fink. Here is the link to the octave-forge on-line documentation:http://octave.sourceforge.net/index/

Contents

Octave Tutorial and Notes:

Create a column vector named: v
```
octave:2> v = [ 0; 1; 2]

v =

0

1
```
Create a row vector named: v
```
octave:1> v = [ 0 1 2 ]

v =

0 1 2
```

Create a 2 x 3 matrix named: m:

octave:3> m = [ 0 1 2; 3 4 5 ]

m =

0 1 2

3 4 5

Create vector containing 11 elements with values between 0 and 10:
```
octave:4> s = [ 0:10]

s =

0 1 2 3 4 5 6 7 8 9 10
```
Create a vector containing 6 sequential even values.
```
octave:5> v = [0:2:10]

v =

0 2 4 6 8 10
```

Create a vector as the composition of two vectors:

octave:6> v = [ 3 4 5 [ 6 7 8 ] ]

v2 =

3 4 5 6 7 8

Create a vector as the composition of two previously created vectors.

octave:8> x0 = [ 6 7 8];

octave:9> x1 = [ 3 2 1 0];

octave:10> x2 = [ x0 x1]

x2 =

6 7 8 3 2 1 0

When creating a vector or matrix using composition the sizes of the component vectors must be compatible. Notice that this example produces an error.

octave:14> x1 = [0; 1; 3] % create a column vector

x1 =

0

1

3

octave:15> x3 = [x0 x1] % compose a row and col vector

error: number of rows must match (3 != 1) near line 15, column 11

error: evaluating assignment expression near line 15, column 4

Use size() and length() to determine the number of rows/columns in a matrix or vector.

octave:15> size(x1)

ans =

3 1

octave:16> size(x0)

ans =

1 3

octave:17> size(m)

ans =

2 3

octave:18> length(x1)

ans = 3

octave:19> length(x0)

ans = 3

Access the values in a matrix or vector by addressing the elements by their index.

 octave:23> disp(v2) % use disp() to display a mtx's values

3 4 5 6 7 8

octave:24> v2(4) % show the 4th element in vector v2

ans = 6

Remember indexing begins with 1. In other words the first element in the vector is at index 1 (not index 0).

octave:25> v2(1)

ans = 3

octave:26> v2(0) % error! index 0 is always invalid

error: invalid vector index = 0

Likewise using an index greater than the length of a vector will produce an error.
```
octave:26> v2(10)

error: invalid vector index = 10
```
Extract a vector of values from a vector. Create a vector by extracting the values from v2 at indexes 2,3 and 4.
```
octave:27> v4 = v2(2:4)

ans =

4 5 6
```
Extract a vector of values from a vector. Create a vector by extracting the values from v2 at indexes 1 and 3. Notice the index argument uses the sequential vector syntax to generate the index values.
```
octave:30> v2(1:2:4)

ans =

3 5
```
Use the 'end' keyword to access all the elements from a start element to the last element in the vector.
```
octave:34> v2(3:end)

ans =

5 6 7 8
```

Extract rows and columns from a matrix.

octave:41> x = [ 0 1 2; 4 5 6; 7 8 9] % create a matrix

x =

0 1 2

4 5 6

7 8 9

octave:42> x(1:2,:) % extract the first two rows

ans =

0 1 2

4 5 6

octave:43> x(:,2:3) % extract columns 2 and 3

ans =

1 2

5 6

8 9

Extract a matrix from a matrix.

octave:44> x(1:2,2:3)

ans =

1 2

5 6

octave:45> x(1:end,2:3)

ans =

1 2

5 6

8 9

Extract all the elements of a matrix as a vector. This is a form of composition. The columns of the vector are concatenated beginning with the first column to form a column vector
```
octave:46> m(:)

ans =

0

3

1

4

2

5
```

Arithmetic operations with vectors and scalars.

octave:53> v0 = [ 1 2 3];

octave:54> v0 + 5 % sum of a vector and a scalar

ans =

6 7 8

octave:55> v0 * 5 % product of a vector and a scalar

ans =

5 10 15

Arithmetic operations on vectors. Note the .* is used to perform an element for element multiply.
```
octave:57> v0 .* [2 4 6]

ans =

2 8 18
```

Here is an example of the error generated when two matrices of the same size are multiplied.

octave:59> v0 * [2 4 6]

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

error: evaluating binary operator `*' near line 59, column 4

To perform an element by element multiply the two vectors must be exactly the same size. Here is an example of the error generated when there is a size mismatch.

octave:59> v0 .* [ 20 30 40 50 ]

error: product: nonconformant arguments (op1 is 1x3, op2 is 1x4)

error: evaluating binary operator `.*' near line 59, column

Matrix transposition: swapping rows and columns of a matrix

octave:60> disp(m)

0 1 2

3 4 5

octave:61> m' % apostrophe operator performs transposition

ans =

0 3

1 4

2 5

Comparison operators (>,<,>=,<=,==,!=) produce matrices of ones and zeros (boolean matrix). In this example the result vector has a one where the value in the vector is greater than 4 and 0 where the value is less than or equal to 4.
```
octave:64> v % display the vector v

v =

0 2 4 6 8 10

octave:65> v > 4 ans =

0 0 0 1 1 1
```
Using the compare operation to select values from a vector. The expression in parenthesis returns [0 0 0 1 1 1] which in turn is multiplied by the vector v. Only those elements which have a corresponding 1 in the boolean vector will occur in the result.
```
octave:69> (v > 4) .* v

ans =

0 0 0 6 8 10
```
Use the find() function to return the indexes whose position evaluates to true.
```
octave:73> disp(v)

0 2 4 6 8 10

octave:74> find(v > 4)

ans =

4 5 6
```
Useful built-in functions: The following built in function will be particularly useful in this class. The versions of wavin() and wavout() linked here will read and write .wav files under OS-X correctly. Download and copy them to your octave folder to use them.
```
size(), length()

help()

max(), min()

zeros(), ones()

rand()

cumsum()

mod()

abs()

pow()

exp()

floor()

plot() and stem()

wavout(), wavin() 
```
Writing octave programs: Below is the simplest possible octave program. It takes two numbers add them together and returns the sum. Remember you must create octave programs using a plain text editor and the scripts must be stored in the folder specified by the LOADPATH variable in your .octaverc file. See the course software installation document for more about this. The name of the file should match the name of the function and use the file name extension .m. This function should therefore be saved as myadder.m.
```
function y = myadder( a, b )

y = a + b;

endfunction
```
Calling an octave program: here is an example of calling the myadder().
```
octave:82> myadder(3,2)

ans = 5
```
Argument passing and return values: When you call an octave program the parameters listed in parenthesis on the first line of the program (a and b in myadder()) are assigned the values in their corresponding position in the function call. In this example a is set to 3 and b is set to 2. The return value (y in myadder()) is assigned the output value which then becomes the value of the function on the command line.
Putting the semicolon after an octave statement prevents the value of the statement from being printed. (The semicolon serves the same purpose on the command line). Try removing it from myadder() to see the result. Note that leaving the semicolon off of statements in a function is a good way of debugging the function because the result of the statements will be printed as they are executed. The printed values will allow you to trace the state of the function as it evaluates.

Multiple values and matrices can be returned from a program. Create a file named summDiff.m with the following code:

function [v0 v1] = sumDiff( arg0, arg1 )

v0 = arg0 + arg1; % create the sum or arg0 and arg1

v1 = arg0 - arg1; % create the difference of arg0 and arg1

endfunction

Accessing multiple return values: Note the syntax used to access the multiple return values. In this case x0 and x1 are assigned values of v0 and v1 from the function.
```
[x0 x1] =sumDiff([1 2 3],[4 5 6])

x0 =

5 7 9

x1 =

-3 -3 -3
```
One reason to write a function is to take advantage of the program flow control that is available there. Two of the most common flow control constructs are if-else-endif and the for loop. This function uses both techniques:
```
function outSig = sGen( srate, frqHz, durSeconds, cosFl )

% sGen( srate, frqHz, durSeconds, cosFl )

% Purpose: calculate a sine signal as a vector

% Input:

% srate = signal sample rate

% frqHz = signal frequency in Hertz

% durSeconds = length of the signal in seconds

% cosFl = 0 to compute a sine 1 to compute a cosine.

% Output: vector containing the a sinusoid of the specified sample rate

% duration and frequency.



% if a cosine was requested then set the initial phase to pi/2

if cosFl != 0

	phaseOffset = pi/2;

else

	% otherwise set initial phase to 0

	phaseOffset = 0;

endif



% calculate the number of samples in the output signal

N = floor(durSeconds * srate);



% initialize the output signal vector to all zeros

outSig = zeros(1,N);



% for each value of n between 0 and N calculate the output sample

for n = 0:N-1

	outSig(n+1) = sin( (2*pi*frqHz*n/srate) + phaseOffset);

endfor



% a computationally faster and more succinct way to compute a sine signal

% outSig = sin( initPhase + (2*pi*frqHz/srate) .* [0:N-1])



endfunction
```
1. The comment block immediately following the first line of the function is given as the help text for the function. At the octave prompt type: help sGen to see how it is used.
2. Notice how the code is indented. The body of the if statement and for loop is indented. This is conventional text programming style and is generally considered to make the code easier to read. The indenting however has no bearing on the functionality of the program.
3. The floor() function is used force the number of samples in the output signal to be an integer. Needless to say it wouldn't make sense to have a fractional count of samples in a signal (What would it mean to have ½ a sample?). Using floor() is a good idea when computing a value which must be an integer – like a sample count or matrix index. Use help floor to find out the algorithm floor() uses to convert fractional number into integers.
4. Other useful flow control statements are: while(), switch(), and elseif. Look these up in the octave documentation to find out how they work.

Contents

Getting Started with Octave Programing

The easiest way to begin writing octave programs is to begin with an existing script and modify it in a series of small increments. This approach is often effective when learning a new programming language. The idea is to start from a simple working function, strip out it's contents, and incrementally replace it with new code. After each small change the program is executed. This way when the program doesn't work as expected or generates an error message it is clear which change produced the error. You can then back out of the change and try an alternative or search through the octave documentation to find a solution.

These directions assume that you have already gone over the tutorial above and feel comfortable with creating, manipulating and displaying octave matrices. If you want more practice you should go through the on-line tutorial also.

Directions for writing a new octave function:

Make a copy of the sumDiff.m script given above.
Rename the file and function to the to the name of the new function you are writing (e.g. myNewFunction) .
From the octave command prompt run the new function. It should behave exactly like the original sumDiff function.
Change the argument list of myNewFunction to the arguments for the new function and delete (or comment out) the body of the function left over from sumDiff.

In the body of myNewFunction print out the values of the new argument list.

function y = myNewFunction( arg0, arg1, arg2 )

	arg0

	arg1

	arg2

	y = 0

endfunction

Run the function again. This time the arguments that you pass should be printed to the console.
```
octave:1> myNewFunction(1,2,[4 5 6])

arg0 = 1

arg1 = 2

arg2 =

4 5 6

y = 10

ans = 10
```
Now begin inserting new code into the body of the function. Don't add more than a couple of lines without running the script again.
- If you get stuck on particular error or aren't sure how to proceed try experimenting at the octave command prompt with simple examples of the octave commands you are using.
- While you are developing signal processing functions start with very low sample rates (e.g. 16 or 32). Low sample rates produce small signal vectors which you can easily verify using the plot() function or by examining the values directly. As a rule listening to the output of an audio signal processing function should only be done after you have verified the process numerically or graphically at a low sample rate.
- Leave the semicolons off the end of each new line of code. This way the result of the calculation is printed to the console when you run it. Once you are sure the calculation is producing a valid value append the semicolon so you are not confused by cluttered output on the console window.
- Many bugs in octave programming are a result of mismatched matrix sizes. If you encounter error messages like this:
```
error: product: nonconformant arguments (op1 is 1x3, op2 is 1x4)

error: evaluating binary operator `.*' near line 59, column 
```
  print out the size of the vector just prior to the operation which generates the message. For example if the the code
```
v0 + v1
```
  is producing the error. Change the code to:
```
size(v0)

size(v1)

v0 + v1
```
  to print out the vector sizes. If the sizes still appear to be correctly matched then examine the documentation for the function or operation you are performing to be sure that you are clear about the data size requirements. Don't forget that the element by element multiply and divide operations on matrices require use of dot version of the multiply and divide operators (e.g. .* or ./ )
Don't believe the line numbers reported in the error message too much. Often the reported line number is just after (or before) the line containing the error. You should treat the line numbers as suggestions regarding where to look for the error rather than as absolute fact. The same is true of the error messages themselves. Often the message will indicate a different problem than the one which actually exists. This shouldn't be too surprising; if octave knew the real problem it would fix it and not bother you about it. This is just another reason to use the incremental development approach. After solving a few bugs you will quickly start to recognize patterns in the way the program responds and reports trouble.

Contents

Numeric Data Types

Clikable LINK.
http://sunsite.univie.ac.at/textbooks/octave/octave_5.html

Go to the first, previous, next, last section, table of contents.

Numeric Data Types

A numeric constant may be a scalar, a vector, or a matrix, and it may contain complex values.

The simplest form of a numeric constant, a scalar, is a single number that can be an integer, a decimal fraction, a number in scientific (exponential) notation, or a complex number. Note that all numeric constants are represented within Octave in double-precision floating point format (complex constants are stored as pairs of double-precision floating point values). Here are some examples of real-valued numeric constants, which all have the same value:

105
1.05e+2
1050e-1

To specify complex constants, you can write an expression of the form

3 + 4i
3.0 + 4.0i
0.3e1 + 40e-1i

all of which are equivalent. The letter `i' in the previous example stands for the pure imaginary constant, defined as

For Octave to recognize a value as the imaginary part of a complex constant, a space must not appear between the number and the `i'. If it does, Octave will print an error message, like this:

octave:13> 3 + 4 i

parse error:

 3 + 4 i
       ^

You may also use `j', `I', or `J' in place of the `i' above. All four forms are equivalent.

Matrices

It is easy to define a matrix of values in Octave. The size of the matrix is determined automatically, so it is not necessary to explicitly state the dimensions. The expression

a = [1, 2; 3, 4]

results in the matrix

Elements of a matrix may be arbitrary expressions, provided that the dimensions all make sense when combining the various pieces. For example, given the above matrix, the expression

[ a, a ]

produces the matrix

but the expression

[ a, 1 ]

produces the error

error: number of rows must match near line 13, column 6

(assuming that this expression was entered as the first thing on line 13, of course).

Inside the square brackets that delimit a matrix expression, Octave looks at the surrounding context to determine whether spaces and newline characters should be converted into element and row separators, or simply ignored, so commands like

[ linspace (1, 2) ]

and

a = [ 1 2
     3 4 ]

will work. However, some possible sources of confusion remain. For example, in the expression

[ 1 - 1 ]

the `-' is treated as a binary operator and the result is the scalar 0, but in the expression

[ 1 -1 ]

the `-' is treated as a unary operator and the result is the vector [ 1, -1 ].

Given a = 1, the expression

[ 1 a' ]

results in the single quote character `'' being treated as a transpose operator and the result is the vector [ 1, 1 ], but the expression

[ 1 a ' ]

produces the error message

error: unterminated string constant

because to not do so would make it impossible to correctly parse the valid expression

[ a 'foo' ]

For clarity, it is probably best to always use commas and semicolons to separate matrix elements and rows. It is possible to enforce this style by setting the built-in variable whitespace_in_literal_matrix to "ignore".

Built-in Variable: whitespace_in_literal_matrix

This variable allows some control over how Octave decides to convert spaces to commas and semicolons in matrix expressions like [m (1)] or

[ 1, 2,
 3, 4 ]

If the value of whitespace_in_literal_matrix is "ignore", Octave will never insert a comma or a semicolon in a literal matrix list. For example, the expression [1 2] will result in an error instead of being treated the same as [1, 2], and the expression

[ 1, 2,
 3, 4 ]

will result in the vector [ 1, 2, 3, 4 ] instead of a matrix.

If the value of whitespace_in_literal_matrix is "traditional", Octave will convert spaces to a comma between identifiers and `('. For example, given the matrix

m = [3 2]

the expression

[m (1)]

will be parsed as

[m, (1)]

and will result in

[3 2 1]

and the expression

[ 1, 2,
 3, 4 ]

will result in a matrix because the newline character is converted to a semicolon (row separator) even though there is a comma at the end of the first line (trailing commas or semicolons are ignored). This is apparently how MATLAB behaves.

Any other value for whitespace_in_literal_matrix results in behavior that is the same as traditional, except that Octave does not convert spaces to a comma between identifiers and `('. For example, the expression

[m (1)]

will produce `3'. This is the way Octave has always behaved.

When you type a matrix or the name of a variable whose value is a matrix, Octave responds by printing the matrix in with neatly aligned rows and columns. If the rows of the matrix are too large to fit on the screen, Octave splits the matrix and displays a header before each section to indicate which columns are being displayed. You can use the following variables to control the format of the output.

Built-in Variable: output_max_field_width: This variable specifies the maximum width of a numeric output field. The default value is 10.

Built-in Variable: output_precision: This variable specifies the minimum number of significant figures to display for numeric output. The default value is 5.

It is possible to achieve a wide range of output styles by using different values of output_precision and output_max_field_width. Reasonable combinations can be set using the format function. See section Basic Input and Output.

Built-in Variable: split_long_rows

For large matrices, Octave may not be able to display all the columns of a given row on one line of your screen. This can result in missing information or output that is nearly impossible to decipher, depending on whether your terminal truncates or wraps long lines.

If the value of split_long_rows is nonzero, Octave will display the matrix in a series of smaller pieces, each of which can fit within the limits of your terminal width. Each set of rows is labeled so that you can easily see which columns are currently being displayed. For example:

octave:13> rand (2,10)
ans =

Columns 1 through 6:

 0.75883  0.93290  0.40064  0.43818  0.94958  0.16467
 0.75697  0.51942  0.40031  0.61784  0.92309  0.40201

Columns 7 through 10:

 0.90174  0.11854  0.72313  0.73326
 0.44672  0.94303  0.56564  0.82150

The default value of split_long_rows is nonzero.

Octave automatically switches to scientific notation when values become very large or very small. This guarantees that you will see several significant figures for every value in a matrix. If you would prefer to see all values in a matrix printed in a fixed point format, you can set the built-in variable fixed_point_format to a nonzero value. But doing so is not recommended, because it can produce output that can easily be misinterpreted.

Built-in Variable: fixed_point_format

If the value of this variable is nonzero, Octave will scale all values in a matrix so that the largest may be written with one leading digit. The scaling factor is printed on the first line of output. For example,

octave:1> logspace (1, 7, 5)'
ans =

 1.0e+07  *

 0.00000
 0.00003
 0.00100
 0.03162
 1.00000

Notice that first value appears to be zero when it is actually 1. For this reason, you should be careful when setting fixed_point_format to a nonzero value.

The default value of fixed_point_format is 0.

Empty Matrices

A matrix may have one or both dimensions zero, and operations on empty matrices are handled as described by Carl de Boor in An Empty Exercise, SIGNUM, Volume 25, pages 2--6, 1990 and C. N. Nett and W. M. Haddad, in A System-Theoretic Appropriate Realization of the Empty Matrix Concept, IEEE Transactions on Automatic Control, Volume 38, Number 5, May 1993.

By default, dimensions of the empty matrix are printed along with the empty matrix symbol, `[]'. The built-in variable print_empty_dimensions controls this behavior.

Built-in Variable: print_empty_dimensions

If the value of print_empty_dimensions is nonzero, the dimensions of empty matrices are printed along with the empty matrix symbol, `[]'. For example, the expression

zeros (3, 0)

will print

ans = [](3x0)

Empty matrices may also be used in assignment statements as a convenient way to delete rows or columns of matrices. See section Assignment Expressions.

Octave will normally issue a warning if it finds an empty matrix in the list of elements that make up another matrix. You can use the variable empty_list_elements_ok to suppress the warning or to treat it as an error.

Built-in Variable: empty_list_elements_ok

This variable controls whether Octave ignores empty matrices in a matrix list.

For example, if the value of empty_list_elements_ok is nonzero, Octave will ignore the empty matrices in the expression

a = [1, [], 3, [], 5]

and the variable a will be assigned the value [ 1, 3, 5 ].

The default value is "warn".

When Octave parses a matrix expression, it examines the elements of the list to determine whether they are all constants. If they are, it replaces the list with a single matrix constant.

Built-in Variable: propagate_empty_matrices: If the value of propagate_empty_matrices is nonzero, functions like inverse and svd will return an empty matrix if they are given one as an argument. The default value is 1.

Ranges

A range is a convenient way to write a row vector with evenly spaced elements. A range expression is defined by the value of the first element in the range, an optional value for the increment between elements, and a maximum value which the elements of the range will not exceed. The base, increment, and limit are separated by colons (the `:' character) and may contain any arithmetic expressions and function calls. If the increment is omitted, it is assumed to be 1. For example, the range

1 : 5

defines the set of values `[ 1, 2, 3, 4, 5 ]', and the range

1 : 3 : 5

defines the set of values `[ 1, 4 ]'.

Although a range constant specifies a row vector, Octave does not convert range constants to vectors unless it is necessary to do so. This allows you to write a constant like `1 : 10000' without using 80,000 bytes of storage on a typical 32-bit workstation.

Note that the upper (or lower, if the increment is negative) bound on the range is not always included in the set of values, and that ranges defined by floating point values can produce surprising results because Octave uses floating point arithmetic to compute the values in the range. If it is important to include the endpoints of a range and the number of elements is known, you should use the linspace function instead (see section Special Utility Matrices).

When Octave parses a range expression, it examines the elements of the expression to determine whether they are all constants. If they are, it replaces the range expression with a single range constant.

Predicates for Numeric Objects

Function File: is_matrix (a): Return 1 if a is a matrix. Otherwise, return 0.

Function File: is_vector (a): Return 1 if a is a vector. Otherwise, return 0.

Function File: is_scalar (a): Return 1 if a is a scalar. Otherwise, return 0.

Function File: is_square (x): If x is a square matrix, then return the dimension of x. Otherwise, return 0.

Function File: is_symmetric (x, tol): If x is symmetric within the tolerance specified by tol, then return the dimension of x. Otherwise, return 0. If tol is omitted, use a tolerance equal to the machine precision.

Go to the first, previous, next, last section, table of contents.

Saturday, December 8, 2007

Using Matlab for Higher Order ODEs and Systems of ODEs

Original document

(Continuation of Using Matlab for First Order ODEs)

Numerical Solution: Converting problems to first order systems
Plotting the solution
Finding numerical values at given t values
Making phase plane plots
Vector fields for autonomous problems: Plotting the vector field
Plotting the vector field together with solution curves
Symbolic Solution: Example with 2nd order problem
Plotting the solution
Finding numerical values at given t values
Example with first order system
Plotting the solution
Finding numerical values at given t values
Making phase plane plots

Numerical solution

Example problem: The angle y of an undamped pendulum with a driving force sin(5 t) satisfies the differential equation

y'' = -sin(y) + sin(5 t)

and the initial conditions

y(0) = 1
y'(0) = 0.

If your problem is of order 2 or higher: rewrite your problem as a first order system. Let y₁=y and y₂=y', this gives the first order system

y₁' = y₂,
y₂' = -sin(y₁) + sin(5 t)

with the initial conditions

y₁(0) = 1
y₂(0) = 0.

Define an inline function g for the right hand side of the first order system: Note that g must be specified as a column vector using [...;...] (not a row vector using [...,...]). In the definition of g, use y(1) for y₁, use y(2) for y₂. The definition of g should have the form

g = inline('[ expression for y₁' ; expression for y₂' ]', 't', 'y');

You have to use inline(...,'t','y'), even if t does not occur in your formula. For our example the first component of g is y₂, the second component is -sin(y₁) + sin(5 t) and we define

g = inline('[ y(2); -sin(y(1))+sin(5*t) ]', 't', 'y')

To plot the numerical solution: To obtain plots of all the components of the solution for t going from t0 to t1 use ode45(g,[t0,t1],[y10;y20]) where y10, y20 are the initial values for y₁, y₂ at the starting point t0. In our example we get a plot for t going from 0 to 20 with the initial values [1;0] using

ode45(g,[0,20],[1;0])

This shows the two functions y₁(t)=y(t) (blue) and y₂(t)=y'(t) (green).

The circles mark the values which were actually computed (the points are chosen by Matlab to optimize accuracy and efficiency). You can obtain a vector ts and a matrix ys with the coordinates of these points using [ts,ys] = ode45(g,[t0,t1],[y10;y20]). You can then plot the solution curves using plot(ts,ys) (this is a way to obtain a plot without the circles). Note that each row of the matrix ys contains 2 entries corresponding to the two components of the solution at that time:

[ts,ys] = ode45(g,[0,20],[1;0]); % find ts, ys, but don't show plot(ts,ys) % make plot of y1 and y2 vs. t [ts,ys] % show table with 3 columns for t, y1, y2

You can obtain the vector of y1 values in the first column of ys by using ys(:,1), therefore plot(ts,ys(:,1)) plots only y₁(t).

To obtain numerical values at specific t values: You can specify a vector tv of t values and use [ts,ys] = ode45(g,tv,[y10;y20]). The first element of the vector tv is the initial t value; the vector tv must have at least 3 elements. E.g., to obtain the solution with the initial values [1;0] at t = 0, 0.5, ..., 10 and display the results as a table with 3 columns t, y₁, y₂, use

[ts,ys]=ode45(g,0:0.5:10,[1;0]); [ts,ys]

To plot trajectories in the phase plane: To see the points with coordinates ( y₁(t), y₂(t) ) in the y₁, y₂ plane for t going from 0 to 20 type

options=odeset('OutputFcn','odephas2'); ode45(g,[0,20],[1;0],options)

This shows the points while they are being computed (the plotting can be stopped with the stop button). To first compute the numerical values and then plot the curve without the circles, type

[ts,ys] = ode45(g,[0,20],[1;0]); % find ts, ys, but don't show plot(ys(:,1),ys(:,2)) % make plot of y2 vs. y1

Vector fields for autonomous systems of two first order ODEs

If the right hand side function g(t, y) does not depend on t, the problem is called autonomous. In this case the behavior of the differential equation can be visualized by plotting the vector g(t, y) at each point y = (y₁,y₂) in the y₁,y₂ plane (the so-called phase plane).

First save the files vectfield.m and vectfieldn.m into the same directory where your m-files are.

To plot the vector field for y1 going from a1 to b1 with a spacing of d1 and y2 going from a2 to b2 with a spacing of d2 use vectfield(g,a1:d1:b1,a2:d2:b2). The command vectfieldn works in the same way, but produces arrows which all have the same length. This makes it easier to see the direction of the vector field.

Example: The undamped pendulum problem without a driving force has the differential equation y'' = -sin(y). This gives the first order system

y₁' = y₂,
y₂' = -sin(y₁)

Here we define

g = inline('[y(2);-sin(y(1))]','t','y')

We can plot the vector field and 10 trajectories with starting points (0,0), (0,0.3), ..., (0,2.7) in the phase plane as follows:

vectfield(g,-2:.5:8,-2.5:.25:2.5)
hold on
for y20=0:0.3:2.7
 [ts,ys] = ode45(g,[0,10],[0;y20]);
 plot(ys(:,1),ys(:,2))
end
hold off

Symbolic solution

Use the dsolve command. Specify all differential equations as strings, using Dy for y'(t), D2y for y''(t) etc. .

For an initial value problem specify the initial conditions in the form 'y(t0)=y0', 'Dy(t0)=y1' etc. . The last argument of dsolve is the name of the independent variable, e.g., 't'.

Example: For the differential equation

y'' = -y + sin(5 t)

use

sol = dsolve('D2y = -y + sin(5*t)','t')

In this case, the answer can be simplified by typing

s = simple(sol)

which gives the general solution -1/24*sin(5*t)+C1*cos(t)+C2*sin(t) with two constants C1, C2.

To solve the ODE with initial conditions y(0) = 1, y'(0) = 0 use

sol = dsolve('D2y = -y + sin(5*t)','y(0)=1','Dy(0)=0','t')

Then s = simple(sol) gives the solution -1/24*sin(5*t)+5/24*sin(t)+cos(t).

To plot the solution curve use ezplot:

ezplot(sol, [0,20])

To obtain numerical values at one or more t values proceed exactly as in the case of a first order ODE.

Example for system of ODEs: For the system

y₁' = y₂,
y₂' = -y₁ + sin(5 t)

with the initial conditions

y₁(0) = 1
y₂(0) = 0

type

sol = dsolve('Dy1=y2','Dy2=-y1+sin(5*t)','y1(0)=1','y2(0)=0','t')

which gives the somewhat strange response

sol = y1: [1x1 sym] y2: [1x1 sym]

This means that the two components of the solution can be accessed as sol.y1 and sol.y2 : Typing

sol.y1

gives -2/3*sin(t)*cos(t)^4+1/2*sin(t)*cos(t)^2+1/6*sin(t)+cos(t). This can be simplified by typing

s1 = simple(sol.y1)

which gives -1/24*sin(5*t)+5/24*sin(t)+cos(t). For the second component y2 of the solution we proceed in the same way: Typing

s2 = simple(sol.y2)

gives -sin(t)-5/24*cos(5*t)+5/24*cos(t).

To plot the solution curves use ezplot:

ezplot(sol.y1, [0,20]) hold on ezplot(sol.y2, [0,20]) hold off

To obtain numerical values use double(subs(sol.y1,'t',tval)) where tval is a number or a vector of numbers. E.g., the following commands generate a table with three columns t, y1, y2 for t=0, 0.5, ..., 10:

tval = (0:0.5:10)'; % column vector with t-values yval = double(subs([sol.y1,sol.y2],'t',tval)); % 2 columns with y1,y2 [tval, yval] % display 3 columns together

To plot the solution in the phase plane use a similar approach with more t values:

tval = (0:0.1:10)'; % column vector with t-values yval = double(subs([sol.y1,sol.y2],'t',tval)); % 2 columns with y1,y2 plot(yval(:,1),yval(:,2)) % plot col.2 of yval vs. col.1 of yval

MATH 246 course page

Tobias von Petersdorff

Using Matlab for First Order ODEs

Original documents

Inline functions
Direction fields
Numerical solution of initial value problems: Plotting the solution
Combining direction field and solution curves
Finding numerical values at given t values
Symbolic solution of ODEs: Finding the general solution
Solving initial value problems
Plotting the solution
Finding numerical values at given t values

Inline Functions

If you want to use a function several times it is convenient to define it as a so-called inline function:

f1 = inline('sin(x)*x','x')

defines the function f1(x)=sin(x)*x. Note that the arguments of inline must be strings (not symbolic expressions). You can then use the function f1 in expressions you type in.

You can also define inline functions of several variables:

g1 = inline('x*y+sin(x)','x','y')

defines the function g1(x,y)=x*y+sin(x) of two variables.

Direction Fields

First download the file dirfield.m and put it in the same directory as your other m-files for the homework.

Define an inline function g of two variables t, y corresponding to the right hand side of the differential equation y'(t) = g(t,y(t)). E.g., for the differential equation y'(t) = t y² define

g = inline('t*y^2','t','y')

You have to use inline(...,'t','y'), even if t or y does not occur in your formula.

To plot the direction field for t going from t0 to t1 with a spacing of dt and y going from y0 to y1 with a spacing of dy use dirfield(g,t0:dt:t1,y0:dy:y1). E.g., for t and y between -2 and 2 with a spacing of 0.2 type

dirfield(g,-2:0.2:2,-2:0.2:2)

Solving an initial value problem numerically

First define the inline function g corresponding to the right hand side of the differential equation y'(t) = g(t,y(t)). E.g., for the differential equation y'(t) = t y² define

g = inline('t*y^2','t','y')

To plot the numerical solution of an initial value problem: For the initial condition y(t0)=y0 you can plot the solution for t going from t0 to t1 using ode45(g,[t0,t1],y0).

Example: To plot the solution of the initial value problem y'(t) = t y², y(-2)=1 in the interval [-2,2] use

ode45(g,[-2,2],1)

The circles mark the values which were actually computed (the points are chosen by Matlab to optimize accuracy and efficiency). You can obtain vectors ts and ys with the coordinates of these points using [ts,ys] = ode45(g,[t0,t1],y0). You can then plot the solution using plot(ts,ys) (this is a way to obtain a plot without the circles).

To combine plots of the direction field and several solution curves use the commands hold on and hold off: After obtaining the first plot type hold on, then all subsequent commands plot in the same window. After the last plot command type hold off.

Example: Plot the direction field and the 13 solution curves with the initial conditions y(-2) = -0.4, -0.2, ..., 1.8, 2:

dirfield(g,-2:0.2:2,-2:0.2:2)
hold on
for y0=-0.4:0.2:2
 [ts,ys] = ode45(g,[-2,2],y0); plot(ts,ys)
end
hold off

To obtain numerical values of the solution at certain t values: You can specify a vector tv of t values and use [ts,ys] = ode45(g,tv,y0). The first element of the vector tv is the initial t value; the vector tv must have at least 3 elements. E.g., to obtain the solution with the initial condition y(-2)=1 at t = -2, -1.5, ..., 1.5, 2 and display the results as a table with two columns, use

[ts,ys]=ode45(g,-2:0.5:2,1); [ts,ys]

To obtain the numerical value of the solution at the final t-value use ys(end) .

Solving a differential equation symbolically

You have to specify the differential equation in a string, using Dy for y'(t) and y for y(t): E.g., for the differential equation y'(t) = t y² type

sol = dsolve('Dy=t*y^2','t')

The last argument 't' is the name of the independent variable. Do not type y(t) instead of y.

If Matlab can't find a solution it will return an empty symbol. If Matlab finds several solutions it returns a vector of solutions.

Sometimes Matlab can't find an explicit solution, but returns the solution in implicit form. E.g., dsolve('Dy=1/(y-exp(y))','t') returns

t-1/2*y^2+exp(y)+C1=0

Unfortunately Matlab cannot handle initial conditions in this case. You can use ezcontour('t-1/2*y^2+exp(y)',[-4 4 -3 3]) to plot several solution curves for t in [-4,4], y in [-3,3]. You can use ezplot('t-1/2*y^2+exp(y)-1',[-4 4 -3 3]) to plot only the curve where t-1/2*y^2+exp(y)=1.

The solution will contain a constant C1. You can substitute values for the constant using subs(sol,'C1',value). E.g., to set C1 to 5 and plot this solution for t=-2 to 2 use

ezplot( subs(sol,'C1',5) , [-2 2] )

To solve an initial value problem additionally specify an initial condition:

sol = dsolve('Dy=t*y^2','y(-2)=1','t')

To plot the solution use ezplot(sol,[t0,t1]). Here is an example for plotting the 13 solution curves with the initial conditions y(-2) = -0.4, -0.2, ..., 1.8, 2:

sol = dsolve('Dy=t*y^2','y(-2)=y0','t')
for y0=-0.4:0.2:2
 ezplot( subs(sol,'y0',y0) , [-2 2])
 hold on
end
hold off
axis tight

To obtain numerical values at one or more t values use subs(sol,'t',tval) and double (or vpa for more digits):

sol = dsolve('Dy=t*y^2','y(-2)=1','t')

This gives a numerical value of the solution at t=0.5:

double( subs(sol,'t',0.5) )

This computes numerical values of the solution at t=-2, -1.5, ..., 2 and displays the result as a table with two columns:

tval = (-2:0.5:2)'; % column vector with t-values yval = double( subs(sol,'t',tval) )% column vector with y-values [tval,yval] % display 2 columns together

(continued in Using Matlab for Higher Order ODEs and Systems of ODEs)

MATH 246 course page

Tobias von Petersdorff

Wednesday, October 10, 2007

Links for Plotting

Annotation - From Loren on the Art of Matlab

Plotting in Matlab

Matlab Commands List

Complex Data (off topic!)

Notes on Matlab Basic and Complex Numbers (pdf file)

Plot help from Matlab (online version)

My Octave/Matlab

Wednesday, May 7, 2008

Print to png file from Octave

Monday, May 5, 2008

Octave Tutorial

Link

Contents:

Octave resources:

Octave Tutorial and Notes:

Getting Started with Octave Programing

Numeric Data Types

Numeric Data Types

Matrices

Empty Matrices

Ranges

Predicates for Numeric Objects

Saturday, December 8, 2007

Using Matlab for Higher Order ODEs and Systems of ODEs

Contents

Numerical solution

Vector fields for autonomous systems of two first order ODEs

Symbolic solution

Using Matlab for First Order ODEs

Contents

Inline Functions

Direction Fields

Solving an initial value problem numerically

Solving a differential equation symbolically

Wednesday, October 10, 2007

Links for Plotting

Tuesday, October 9, 2007

Basic Matrix, collected

17.1 Basic Matrix Functions, from GNU Octave Manual by John W. Eaton

Octave Examples for Solving Linear Algebra Equations, from Hanson's website.

Note: very useful for the Intro to Numerical Analysis course.

Basic Octave

General

Blog Archive