Cookbook: Difference between revisions

From Octave
Jump to navigation Jump to search
(→‎Programs, Libraries, and Packages: add recipe to check if package is missing)
 
(10 intermediate revisions by 7 users not shown)
Line 1: Line 1:
An Octave cookbook. Each entry should go in a separate section and have the following subsection: problem, solution, discussion and maybe a see also.
An Octave cookbook. Each entry should go in a separate section and have the following subsection: problem, solution, discussion and maybe a see also.
__FORCETOC__


== Programs, Libraries, and Packages ==
== Programs, Libraries, and Packages ==
Line 15: Line 17:
This is a list of possible tests to check for features:
This is a list of possible tests to check for features:


    ## support for 64 bit indexing
<syntaxhighlight lang="Octave">
    sizemax () > intmax ("int32")
## support for 64 bit indexing
sizemax () > intmax ("int32")


    ## built with support for java
## built with support for java
    usejava ("jvm")
usejava ("jvm")


    ## Image IO with support for tif files
## Image IO with support for tif files
    any (cellfun (@(x) ismember ("tif", x), {imformats.ext}))
any (cellfun (@(x) ismember ("tif", x), {imformats.ext}))
    ## Image IO with support for png files
## Image IO with support for png files
    any (cellfun (@(x) ismember ("png", x), {imformats.ext}))
any (cellfun (@(x) ismember ("png", x), {imformats.ext}))
</syntaxhighlight>


=== Find if a package is installed ===
=== Find if a package is installed ===
Line 35: Line 39:
==== Solution ====
==== Solution ====


Use {{codeline|pkg ("describe", pkg-name)}} like so:
Use {{codeline|pkg ("list", pkg-name)}} like so:


if (! isempty (pkg ("describe", "foo")))
<syntaxhighlight lang="Octave">
  ## use functions from package foo, the prefered way
if (! isempty (pkg ("list", "foo")))
elseif (! isempty (pkg ("describe", "bar")))
  ## use functions from package foo, the preferred way
  ## use functions from package bar, not so optimal
elseif (! isempty (pkg ("list", "bar")))
else
  ## use functions from package bar, not so optimal
  ## default case
else
endif
  ## default case
endif
</syntaxhighlight>


==== Discussion ====
==== Discussion ====
Line 56: Line 62:
then you might as well not catch it in the first place.
then you might as well not catch it in the first place.


## This contraption doesn't add anything.  If 'pkg load' fails, it
<syntaxhighlight lang="Octave">
## will already give an error message to the user.
## This contraption doesn't add anything.  If 'pkg load' fails, it
try
## will already give an error message to the user.
  pkg load foo;
try
catch
  pkg load foo;
  error ("failed to load foo: %s", lasterr ());
catch
end_try_catch
  error ("failed to load foo: %s", lasterr ());
end_try_catch


## Again, doesn't add anything.  The failure of 'pkg load' is enough
## Again, doesn't add anything.  The failure of 'pkg load' is enough
if (isempty (pkg ("describe", "foo")))
if (isempty (pkg ("list", "foo")))
  error ("program: package foo is not installed");
  error ("program: package foo is not installed");
endif
endif
</syntaxhighlight>


Beware that an installed package is not always a guarantee that a function
Beware that an installed package is not always a guarantee that a function
will be available.  Some packages may disable functions at build time, or
will be available.  Some packages may disable functions at build time, or
specific functions may have specific runtime requirements.
specific functions may have specific runtime requirements.
=== List all function in Octave ===
Use the following script (filename <code>list_func.m</code>)
<syntaxhighlight lang="Octave">
## List of all builtin (C++) functions and m-file functions
funcs = vertcat (__builtins__ (), __list_functions__ ());
## Write list to file
fid = fopen ("all_funcs.tmp", "w");
if (fid == -1)
  error ("Unable to open temporary file all_funcs.tmp.  Aborting...\n");
endif
fprintf (fid, "%s\n", funcs{:});
fclose (fid);
</syntaxhighlight>
And execute with
    run-octave -f list_func.m


== Structures ==
== Structures ==
=== Retrieve a field value from all entries in a struct array ===
=== Retrieve a field value from all entries in a struct array ===
==== Problem ====
==== Problem ====
You have a struct array with multiple fields, and you want to access the value from a specific field from all elements. For example, you want to return the age from all patients in the following case:
You have a struct array with multiple fields, and you want to access the value from a specific field from all elements. For example, you want to return the age from all patients in the following case:


  samples = struct ("patient", {"Bob", "Kevin", "Bob" , "Andrew"},
<syntaxhighlight lang="Octave">
                    "age",    { 45  ,  52    ,  45  ,  23    },
samples = struct ("patient", {"Bob", "Kevin", "Bob" , "Andrew"},
                    "protein", {"H2B", "CDK2" , "CDK2", "Tip60" },
                  "age",    { 45  ,  52    ,  45  ,  23    },
                    "tube"  , { 3  ,  5    ,  2    ,  18    }
                  "protein", {"H2B", "CDK2" , "CDK2", "Tip60" },
                    );
                  "tube"  , { 3  ,  5    ,  2    ,  18    }
                  );
</syntaxhighlight>


==== Solution ====
==== Solution ====
Indexing the struct returns a comma separated list so use them to create a matrix.
Indexing the struct returns a comma separated list so use them to create a matrix.


  [samples(:).age]
<syntaxhighlight lang="Octave">
[samples(:).age]
</syntaxhighlight>


This however does not keep the original structure of the data, instead returning all values in a single column. To fix this, use {{Codeline|reshape()}}.
This however does not keep the original structure of the data, instead returning all values in a single column. To fix this, use {{Codeline|reshape()}}.


  reshape ([samples(:).age], size (samples))
<syntaxhighlight lang="Octave">
reshape ([samples(:).age], size (samples))
</syntaxhighlight>


==== Discussion ====
==== Discussion ====
Returning all values in a comma separated lists allows you to make anything out of them. If numbers are expected, create a matrix by enclosing them in square brackets. But if strings are to be expected, a cell array can also be easily generated with curly brackets


  {samples(:).patient}
Returning all values in a comma separated lists allows you to make anything out of them.
If numbers are expected, create a matrix by enclosing them in square brackets.
But if strings are to be expected, a cell array can also be easily generated with curly brackets
 
<syntaxhighlight lang="Octave">
{samples(:).patient}
</syntaxhighlight>


You are also not limited to return all elements, you may use logical indexing from other fields to get values from the others:
You are also not limited to return all elements, you may use logical indexing from other fields to get values from the others:


  [samples([samples(:).age] > 34).tube]                ## return tube numbers from all samples from patients older than 34
<syntaxhighlight lang="Octave">
  [samples(strcmp({samples(:).protein}, "CDK2")).tube]  ## return all tube numbers for protein CDK2
[samples([samples(:).age] > 34).tube]                ## return tube numbers from all samples from patients older than 34
[samples(strcmp({samples(:).protein}, "CDK2")).tube]  ## return all tube numbers for protein CDK2
</syntaxhighlight>


== Array manipulation ==
== Array manipulation ==
Line 109: Line 152:
==== Problem ====
==== Problem ====


For an array {{Codeline|A}} with arbitrary number of dimensions, select, for example, the first column. This would be {{Codeline|A(:, 1)}} if {{Codeline|A}} was 2-D, {{Codeline|A(:, 1, :)}} if {{Codeline|A}} was 3-D, and so on.
For an array {{Codeline|A}} with arbitrary number of dimensions, select, for example, the first column.
This would be {{Codeline|A(:, 1)}} if {{Codeline|A}} was 2-D, {{Codeline|A(:, 1, :)}} if {{Codeline|A}} was 3-D, and so on.


==== Solution ====
==== Solution ====


One possibility is to use {{Codeline|subsref}} with the input {{Codeline|idx}} created dynamically with {{Codeline|repelems}} to have the right number of dimensions. This can be written as a function:
One possibility is to use {{manual|subsref}} with the input {{Codeline|idx}} created dynamically with {{manual|repelems}} to have the right number of dimensions.
 
This can be written as a function:


{{Code||<syntaxhighlight lang="octave" style="font-size:13px">
<syntaxhighlight lang="Octave">
function [B]= array_slice (A,k,d)
function [B]= array_slice (A,k,d)
#return the k-th slice (row, column...) of A, with d specifying the dimension to slice on
#return the k-th slice (row, column...) of A, with d specifying the dimension to slice on
Line 135: Line 179:
%! A=rand(2, 3, 4, 5, 6);
%! A=rand(2, 3, 4, 5, 6);
%!assert (array_slice (A,2,3), A(:, :, 2, :, :))
%!assert (array_slice (A,2,3), A(:, :, 2, :, :))
</syntaxhighlight>}}
</syntaxhighlight>


To remove the singleton dimension {{Codeline|d}} from the result {{Codeline|B}}, use  
To remove the singleton dimension {{Codeline|d}} from the result {{Codeline|B}}, use  
{{Code||<syntaxhighlight lang="octave" style="font-size:13px">
 
<syntaxhighlight lang="Octave">
B = reshape(B, [size(B)([1:d-1 d+1:end])]);
B = reshape(B, [size(B)([1:d-1 d+1:end])]);
</syntaxhighlight>}}
</syntaxhighlight>


== Input/output ==
== Input/output ==
=== Display matched elements from different arrays ===
=== Display matched elements from different arrays ===
==== Problem ====
==== Problem ====
You have two, or more, arrays with paired elements and want to print out a string about them. For example:
You have two, or more, arrays with paired elements and want to print out a string about them. For example:


    keys  = {"human",  "mouse", "chicken"};
<syntaxhighlight lang="Octave">
    values = [ 64        72      70      ];
keys  = {"human",  "mouse", "chicken"};
values = [ 64        72      70      ];
</syntaxhighlight>


and you want to display:
and you want to display:


    Calculated human genome GC content is 64%
Calculated human genome GC content is 64%
    Calculated mouse genome GC content is 72%
Calculated mouse genome GC content is 72%
    Calculated chicken genome GC content is 70%
Calculated chicken genome GC content is 70%


==== Solution ====
==== Solution ====
Make a two rows cell array, with each paired data in a column and supply a cs-list to printf


    values = num2cell (values);
Make a two rows cell array, with each paired data in a column and supply a cs-list to {{manual|printf}}
    new    = {keys{:}; values{:}};
 
    printf ("Calculated %s genome GC content is %i%%\n", new{:})
<syntaxhighlight lang="Octave">
values = num2cell (values);
new    = {keys{:}; values{:}};
printf ("Calculated %s genome GC content is %i%%\n", new{:})
</syntaxhighlight>


or in a single line:
or in a single line:
    printf ("Calculated %s genome GC content is %i%%\n", {keys{:}; num2cell(values){:}}{:})
 
   
<syntaxhighlight lang="Octave">
printf ("Calculated %s genome GC content is %i%%\n", {keys{:}; num2cell(values){:}}{:})
</syntaxhighlight>   


==== Discussion ====
==== Discussion ====
{{Codeline|printf}} and family do not accept cell arrays as values. However, they keep repeating the template given as long as it has enough arguments to keep going. As such, the trick is on supplying a cs-list of elements which can be done by using a cell array and index it with {{Codeline|<nowiki>{}</nowiki>}}.


Since values are stored in column-major order, paired values need to be on the same column. A new row of data can then be added later with {{Codeline|new(end+1,:) <nowiki>= {"Andrew", "Bob", "Kevin"}</nowiki>}}. Note that normal brackets are now being used for indexing.
{{manual|printf}} and family do not accept cell arrays as values.
However, they keep repeating the template given as long as it has enough arguments to keep going.
As such, the trick is on supplying a cs-list of elements which can be done by using a cell array and index it with <code>{}</code>.
 
Since values are stored in column-major order, paired values need to be on the same column.
A new row of data can then be added later with
 
<syntaxhighlight lang="Octave">
new(end+1,:) = {"Andrew", "Bob", "Kevin"};
</syntaxhighlight>
 
Note that normal brackets are now being used for indexing.


=== Swap values ===
=== Swap values ===
If you want to exchange the value between two variables without creating a dummy one, you can simply do:
If you want to exchange the value between two variables without creating a dummy one, you can simply do:


{{Code|Swap values without dummy variable|<syntaxhighlight lang="octave" style="font-size:13px">
<syntaxhighlight lang="Octave">
[b,a] = deal (a,b);
[b,a] = deal (a,b);
</syntaxhighlight>}}
</syntaxhighlight>


=== Collect all output arguments of a function ===
=== Collect all output arguments of a function ===
If you have a function that returns several values, e.g.  
If you have a function that returns several values, e.g.  


{{Code||<syntaxhighlight lang="octave" style="font-size:13px">
<syntaxhighlight lang="Octave">
function [a b c]= myfunc ()
function [a b c]= myfunc ()
   [a,b,c] = deal (1,2,3);  
   [a,b,c] = deal (1,2,3);  
endfunction
endfunction
</syntaxhighlight>}}
</syntaxhighlight>


and you want to collect them all into a single cell (similarly to Python's zip() function) you can do:
and you want to collect them all into a single cell (similarly to Python's zip() function) you can do:


{{Code|Collect multiple output arguments|<syntaxhighlight lang="octave" style="font-size:13px">
<syntaxhighlight lang="Octave">
outargs = nthargout (1:3, @myfunc)
outargs = nthargout (1:3, @myfunc)
</syntaxhighlight>}}
</syntaxhighlight>


=== Create a text table with fprintf===
=== Create a text table with fprintf===
(a.k.a. A funny formatting trick with fprintf found by chance)


Imagine that you want to create a text table with fprintf with 2 columns of 15 characters width and both right justified. How to do this thing?
Imagine that you want to create a text table with {{manual|fprintf}} with 2 columns of 15 characters width and both right justified.
How to do this thing?


That's easy:
That's easy:


If the variable Text is a cell array of strings (of length <15) with two columns and a certain number of rows, simply type for the kth row of Text
If the variable Text is a cell array of strings (of length < 15) with two columns and a certain number of rows, simply type for the k-th row of Text
{{Code||<syntaxhighlight lang="octave" style="font-size:13px">
 
fprintf('%15.15s | %15.15s\n', Text{k,1}, Text{k,2});
<syntaxhighlight lang="Octave">
</syntaxhighlight>}}
fprintf('%15.15s | %15.15s\n', Text{k,1}, Text{k,2});
The syntax '%<n>.<m>s' allocates '<n>' places to write chars and display the '<m>' first characters of the string to display.
</syntaxhighlight>
 
The syntax <code>%<n>.<m>s</code> allocates <code>n</code> places to write chars and display the <code>m</code> first characters of the string to display.


Example:
Example:
{{Code|Example create a text table with fprintf|<syntaxhighlight lang="octave" style="font-size:13px">
 
octave:1> Text={'Hello','World'};
<syntaxhighlight lang="Octave">
octave:2> fprintf('%15.15s | %15.15s\n', Text{1,1}, Text{1,2})
Text = {"Hello", "World"};
fprintf('%15.15s | %15.15s\n', Text{1,1}, Text{1,2})
</syntaxhighlight>
 
           Hello |          World
           Hello |          World
</syntaxhighlight>}}


===Load comma separated values (*.csv) files===
===Load comma separated values (*.csv) files===


{{Code|Load comma separated values files|<syntaxhighlight lang="octave" style="font-size:13px">
# Using {{manual|textread}} gets you a one column cell array.  The original size can be restored by using the {{manual|reshape}} function. <syntaxhighlight lang="Octave">
A=textread("file.csv", "%d", "delimiter", ",");
A = textread("file.csv", "%d", "delimiter", ",");
B=textread("file.csv", "%s", "delimiter", ",");
B = textread("file.csv", "%s", "delimiter", ",");
inds = isnan(A);
inds = isnan(A);
B(!inds) = num2cell(A(!inds))
B(! inds) = num2cell (A(! inds))
</syntaxhighlight>}}
</syntaxhighlight>
# Another option is to use the function {{manual|csvread}}.  However, this function can't handle non-numerical data.
# The probably best option is to use the function [https://octave.sourceforge.io/io/function/csv2cell.html <code>csv2cell</code>] from the [[IO package]].  This function can read mixed-type (numerical and text) .csv files, allows to specify other field separators than a comma and other text protection characters (default: <code>"</code> double quote) and can skip header lines.  If you have the [[IO package]] installed and loaded, type <code>help csv2cell</code> at the Octave prompt for more info.


This gets you a 1 column cell array. You can reshape it to the original size by using the <code>reshape</code> function
===Load XML files===


The next version of octave (3.6) implements the <code>CollectOutput</code> switch as seen in example 8 here:                    http://www.mathworks.com/help/techdoc/ref/textscan.html
Reading XML in octave can be achieved using the java library [https://xerces.apache.org/ Apache Xerces].  


Another option is to use the function <code>csvread</code>, however this function can't handle non-numerical data.
It seems that the Matlab's <code>xmlread</code> is just a thin wrapper around the Apache Xerces library.
One should note however, that Java functions have the working directory set to the working directory when octave starts and the working directory is not modified by a {{manual|cd}} in octave.
Matlab has the same behavior, as Java does not provide a way to change the current working directory (http://bugs.java.com/bugdatabase/view_bug.do?bug_id=4045688).
To avoid any issues, it is thus better to use the absolute path to the XML file.


The probably best option is to use the function csv2cell() from the io package. This function can read mixed-type (numerical and text) .csv files, allows to specify other field separators than a comma and other text protection characters (default: " double quote) and can skip headerlines.  If you have the io package installed and loaded, type "help csv2cell" at the Octave prompt for more info.
You need the jar files {{Path|xercesImpl.jar}} and {{Path|xml-apis.jar}} from e.g. https://xerces.apache.org/mirrors.cgi#binary (check for the latest version).
Use {{manual|javaaddpath}} to include these files:


===Load XML files===
<syntaxhighlight lang="Octave">
javaaddpath ("/path/to/xerces-2_11_0/xercesImpl.jar");
javaaddpath ("/path/to/xerces-2_11_0/xml-apis.jar");
</syntaxhighlight>


Reading XML in octave can be achieved using the java library [http://xerces.apache.org/ Xerces] (from apache).
Example:
 
It seems that the matlab's xmlread is just a thin wrapper around the Xerces library. One should note however, that Java functions have the working directory set to the working directory when octave starts and the working directory is not modified by a cd in octave. Matlab has the same behavior, as Java does not provide a way to change the current working directory (http://bugs.java.com/bugdatabase/view_bug.do?bug_id=4045688). To avoid any issues, it is thus better to use the absolute path to the XML file.


You need the jar files xercesImpl.jar and xml-apis.jar from e.g. https://www.apache.org/dist/xerces/j/Xerces-J-bin.2.11.0.tar.gz (check for the latest version).
<syntaxhighlight lang="Octave">
Use javaaddpath to include these files:
filename = "sample.xml";


{{Code|Define java path|<syntaxhighlight lang="octave" style="font-size:1.1em">
## These three lines are equivalent to xDoc = xmlread(filename) in Matlab
  javaaddpath('/path/to/xerces-2_11_0/xercesImpl.jar');
parser = javaObject("org.apache.xerces.parsers.DOMParser");
  javaaddpath('/path/to/xerces-2_11_0/xml-apis.jar');
parser.parse(filename);  
</syntaxhighlight>}}
xDoc = parser.getDocument();


A sample script:
elem = xDoc.getElementsByTagName("data").item(0); ## get first data element
data = elem.getFirstChild.getTextContent();      ## get text from child
att  = elem.getAttribute("att");                  ## get attribute named att
</syntaxhighlight>


{{Code|Load XML file|<syntaxhighlight lang="octave" style="font-size:1.1em">
{{File|sample.xml|
filename = 'sample.xml';
<syntaxhighlight lang="xml">
<root>
% These 3 lines are equivalent to xDoc = xmlread(filename) in matlab
  <data att="1">hello</data>
parser = javaObject('org.apache.xerces.parsers.DOMParser');
</root>
parser.parse(filename);
xDoc = parser.getDocument;
% get first data element
elem = xDoc.getElementsByTagName('data').item(0);
% get text from child
data = elem.getFirstChild.getTextContent
% get attribute named att
att = elem.getAttribute('att')
</syntaxhighlight>}}
 
The file <tt>sample.xml</tt>:
 
{{Code|Sample XML file|<syntaxhighlight lang="xml" style="font-size:1.1em">
<root>
<data att="1">hello</data>
</root>
</syntaxhighlight>}}
</syntaxhighlight>}}


Line 275: Line 339:
For example, to plot data using a string variable as a legend:
For example, to plot data using a string variable as a legend:


Option 1 (simplest):
# Static string as legend (simplest): <syntaxhighlight lang="Octave">
{{Code|Using variable strings in commands. op1|<syntaxhighlight lang="octave" style="font-size:13px">
x = linspace (-1, 3, 100);
y = sin (x);
legend = "-1;My data;";
legend = "-1;My data;";
plot(x, y, legend);
plot (x, y, legend);
</syntaxhighlight>}}
</syntaxhighlight>
 
# Variable string as legend (moderate): <syntaxhighlight lang="Octave">
Option 2 (to insert variables):
x = linspace (-1, 3, 100);
{{Code|Using variable strings in commands. op2|<syntaxhighlight lang="octave" style="font-size:13px">
y = sin (x);
plot(x, y, sprintf("-1;%s;", dataName));
dataName = "My data";
</syntaxhighlight>}}
plot (x, y, sprintf("-1;%s;", dataName));
 
</syntaxhighlight>
Option 3 (not as neat):
# Variable string as legend using {{manual|eval}} (not as neat): <syntaxhighlight lang="Octave">
{{Code|Using variable strings in commands. op3|<syntaxhighlight lang="octave" style="font-size:13px">
legend = "My data";
legend = 'my legend';
plot_command = ["plot (x, y, '-1;", legend, ";')"];
plot_command = ['plot(x,y,\';',legend,';\')'];
eval (plot_command);
eval(plot_command);
</syntaxhighlight>
</syntaxhighlight>}}


These same tricks are useful for reading and writing data files with unique names, etc.
These same tricks are useful for reading and writing data files with unique names, etc.


== Combinatorics ==
== Combinatorics ==
=== Combinations with string characters ===
=== Combinations with string characters ===
==== Problem ====
==== Problem ====
You want to get all combinations of different letters but {{codeline|nchoosek}} only accepts numeric input.
 
You want to get all combinations of different letters but {{manual|nchoosek}} only accepts numeric input.


==== Solution ====
==== Solution ====
Convert your string to numbers and then back to characters.
Convert your string to numbers and then back to characters.


{{Code||<syntaxhighlight lang="octave">
<syntaxhighlight lang="Octave">
char (nchoosek (uint8 (string), n)
string = "Hello";
</syntaxhighlight>}}
n = 4;
char (nchoosek (uint8 (string), n))
</syntaxhighlight>
 
==== Discussion ====
==== Discussion ====
A string in Octave is just a character matrix and can easily be converted to numeric form back and forth. Each character has an associated number (the {{codeline|asci}} function of the {{forge|miscellaneous}} package displays a nicely formatted conversion table).
 
A string in Octave is just a character matrix and can easily be converted to numeric form back and forth.
Each character has an associated number (the {{codeline|asci}} function of the {{forge|miscellaneous}} package displays a nicely formatted conversion table).


=== Permutations with repetition ===
=== Permutations with repetition ===
==== Problem ====
==== Problem ====
You want to generate all possible permutations of a vector with repetition.
You want to generate all possible permutations of a vector with repetition.


==== Solution ====
==== Solution ====
Use {{codeline|ndgrid}}


{{Code||<syntaxhighlight lang="octave">
Use {{manual|ndgrid}}:
[x y z] = ndgrid ([1 2 3 4 5]);
 
[x(:) y(:) z(:)]
<syntaxhighlight lang="Octave">
</syntaxhighlight>}}
[x, y, z] = ndgrid ([1, 2, 3, 4, 5]);
[x(:), y(:), z(:)]
</syntaxhighlight>


==== Discussion ====
==== Discussion ====
It is possible to expand the code above and make it work for any length of permutations.
It is possible to expand the code above and make it work for any length of permutations.


{{Code||<syntaxhighlight lang="octave">
<syntaxhighlight lang="Octave">
cart  = nthargout ([1:n], @ndgrid, vector);
cart  = nthargout ([1:n], @ndgrid, vector);
combs = cell2mat (cellfun (@(c) c(:), cart, "UniformOutput", false));
combs = cell2mat (cellfun (@(c) c(:), cart, "UniformOutput", false));
</syntaxhighlight>}}
</syntaxhighlight>


== Mathematics ==
== Mathematics ==
=== Test if a number is a integer ===
 
There are several methods to do this. The simplest method is probably {{Codeline|<nowiki>fix (x) == x</nowiki>}}
=== Test if a number is an integer ===
 
The simplest method is probably {{manual|fix}}:
 
<syntaxhighlight lang="Octave">
fix (x) == x
</syntaxhighlight>


=== Find if a number is even/odd ===
=== Find if a number is even/odd ===
==== Problem ====
==== Problem ====
You have a number, or an array or matrix of them, and want to know if any of them is an odd or even number, i.e., their parity.
You have a number, or an array or matrix of them, and want to know if any of them is an odd or even number, i.e., their parity.
==== Solution ====
==== Solution ====
Check the remainder of a division by two. If the remainder is zero, the number is odd.


  mod (value, 2) ## 1 if odd, zero if even
Check the remainder of a division by two. If the remainder is zero, the number is even.


Since {{Codeline|mod()}} acceps a matrix, the following can be done:
<syntaxhighlight lang="Octave">
mod (value, 2) ## 1 if odd, zero if even
</syntaxhighlight>


  any  (mod (values, 2)) ## true if at least one number in values is even
Since {{manual|mod}} accepts a matrix, the following can be done:
  all  (mod (values, 2)) ## true if all numbers in values are odd
 
 
<syntaxhighlight lang="Octave">
  any (!logical (mod (values, 2))) ## true if at least one number in values is even
any  (mod (values, 2)) ## true if at least one number in values is even
  all (!logical (mod (values, 2))) ## true if all numbers in values are even
all  (mod (values, 2)) ## true if all numbers in values are odd
 
any (!logical (mod (values, 2))) ## true if at least one number in values is even
all (!logical (mod (values, 2))) ## true if all numbers in values are even
</syntaxhighlight>


==== Discussion ====
==== Discussion ====
Since we are checking for the remainder of a division, the first choice would be to use {{Codeline|rem()}}. However, in the case of negative numbers {{Codeline|mod()}} will still return a positive number making it easier for comparisons. Another alternative is to use {{Codeline|bitand (X, 1)}} or {{Codeline|bitget (X, 1)}} but those are a bit slower.


Note that this solution applies to integers only. Non-integers such as 1/2 or 4.201 are neither even nor odd. If the source of the numbers are unknown, such as user input, some sort of checking should be applied for NaN, Inf, or non-integer values.
Since we are checking for the remainder of a division, the first choice would be to use {{manual|rem}}.
However, in the case of negative numbers {{manual|mod}} will still return a positive number making it easier for comparisons.
Another alternative is to use {{Codeline|bitand (X, 1)}} or {{Codeline|bitget (X, 1)}} but those are a bit slower.


==== See also ====
Note that this solution applies to integers only.
Find if a number is an integer.
Non-integers such as <code>0.5</code> or <code>4.201</code> are neither even nor odd.
If the source of the numbers are unknown, such as user input, some sort of checking should be applied for <code>NaN</code>, <code>Inf</code>, or non-integer values.


=== Parametrized Functions ===
=== Parametrized Functions ===
==== Problem  ====
==== Problem  ====


Line 365: Line 459:


==== Solution ====
==== Solution ====
We could solve the problem with the following code:


{{Code|Solve spring equation for different values of the spring constant|<syntaxhighlight lang="octave" style="font-size:13px">
We could solve the problem with the following code to solve the spring equation for different values of the spring constant:
 
<syntaxhighlight lang="Octave">
t = linspace (0, 10, 100);
t = linspace (0, 10, 100);
function sprime = spring (s, t, k)
function sprime = spring (s, t, k)
Line 381: Line 476:
plot (t, x1, t, x2)
plot (t, x1, t, x2)
legend ('x1', 'x2')
legend ('x1', 'x2')
</syntaxhighlight>}}
</syntaxhighlight>


[[File:solparfun.png|400px]]
[[File:solparfun.png|400px]]
Line 390: Line 485:


The [http://www.gnu.org/software/octave/doc/interpreter/Anonymous-Functions.html#Anonymous-Functions anonymous function]
The [http://www.gnu.org/software/octave/doc/interpreter/Anonymous-Functions.html#Anonymous-Functions anonymous function]
{{Code||<syntaxhighlight lang="octave" style="font-size:13px">
  @(x, t) sprime (x, t, k)
</syntaxhighlight>}}


is a function of only <math>x, t</math> where the parameter <math>k</math> is 'frozen' to the value it has at the moment in the current scope.
<syntaxhighlight lang="Octave">
@(x, t) sprime (x, t, k)
</syntaxhighlight>
 
is a function of only <math>x, t</math> where the parameter <math>k</math> is "frozen" to the value it has at the moment in the current scope.


=== Distance between points ===
=== Distance between points ===
==== Problem  ====
==== Problem  ====
Given a set of points in space we want to calculate the distance between all of them. Each point is described by its components <math> (x_i,y_i,\ldots)</math>. Asusme that the points are saved in a matrix '''<tt>P</tt>''' with '''<tt>N</tt>''' rows (one for each point) and '''<tt>D</tt>''' columns, one for each component.
 
Given a set of points in space we want to calculate the distance between all of them.
Each point is described by its components <math> (x_i,y_i,\ldots)</math>.
Assume that the points are saved in a matrix <code>P</code> with <code>m</code> rows (one for each point) and <code>n</code> columns, one for each component.
 
==== Solution ====
==== Solution ====
One way of proceeding is to use the broadcast properties of operators in GNU Octave. The square distance between the points can be calculated with the code
 
<!-- {{SyntaxHighlight| -->
One way of proceeding is to use the broadcast properties of operators in GNU Octave.
{{Code|Calculate square distance between points|<syntaxhighlight lang="octave" style="font-size:13px">
The square distance between the points can be calculated with the code:
[N, dim] = size (P);
 
Dsq    = zeros (N);  
<syntaxhighlight lang="Octave">
for i = 1:dim
[m, n] = size (P);
Dsq    = zeros (m);  
for i = 1:n
   Dsq += (P(:,i) - P(:,i)').^2;
   Dsq += (P(:,i) - P(:,i)').^2;
endfor
endfor
</syntaxhighlight>}}
</syntaxhighlight>
 
This matrix is symmetric with zero diagonal.
This matrix is symmetric with zero diagonal.


Similarly the vectors pointing from one point to the another is
Similarly the vectors pointing from one point to the another is
<!-- {{SyntaxHighlight| -->
 
{{Code|Calculate radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px">
<syntaxhighlight lang="Octave">
R     = zeros (N,N,dim);  
R = zeros (m, m, n);  
for i = 1:dim
for i = 1:n
   R(:,:,i) = P(:,i) - P(:,i)';
   R(:,:,i) = P(:,i) - P(:,i)';
endfor
endfor
</syntaxhighlight>}}
</syntaxhighlight>
 
The relation between <code>Dsq</code> and <code>R</code> is


The relation between <tt>Dsq</tt> and <tt>R</tt> is
<syntaxhighlight lang="Octave">
<!-- {{SyntaxHighlight| -->
Dsq = sumsq (R, 3);  
{{Code||<syntaxhighlight lang="octave" style="font-size:13px">
</syntaxhighlight>
Dsq = sumsq (R,3);  
</syntaxhighlight>}}


==== Discussion ====
==== Discussion ====
The calculation can be implemented using functions like <tt>cellfun</tt> and avoid the loop over components of the points. However in most cases we will have more points than components and the improvement, if any, will be minimal.


Another observation is that the matrix Dsq is symmetric and we could store only the lower or upper triangular part. To use this optimization in a practical way check the help of the functions <tt>vech</tt> and <tt>unvech</tt> (this one is in the Forge package ''general''). Two functions that haven't seen the light yet are <tt>sub2ind_tril</tt> and <tt>ind2sub_tril</tt> (currently private functions in the [[Mechanics_package | Forge package mechanics]]) that are useful to index the elements of a vector constructed with the function <tt>vech</tt>. Each page (the third index) of the multidimensional array <tt>R</tt> is an anti-symmetric matrix and we could also save some memory by keeping only one of the triangular submatrices.
The calculation can be implemented using functions like {{manual|cellfun}} and avoid the loop over components of the points.
However in most cases we will have more points than components and the improvement, if any, will be minimal.
 
Another observation is that the matrix <code>Dsq</code> is symmetric and we could store only the lower or upper triangular part.
To use this optimization in a practical way check the help of the functions <code>vech</code> and <code>unvech</code> (this one is in the Forge package {{Forge|general}}).
Two functions that haven't seen the light yet are <code>sub2ind_tril</code> and <code>ind2sub_tril</code> (currently private functions in the [[Mechanics_package | Forge package mechanics]]) that are useful to index the elements of a vector constructed with the function <code>vech</code>.
Each page (the third index) of the multidimensional array <code>R</code> is an anti-symmetric matrix and we could also save some memory by keeping only one of the triangular submatrices.


Check the [[Geometry package]] for many more distance functions (points, lines, polygons, etc.).
Check the [[Geometry package]] for many more distance functions (points, lines, polygons, etc.).


=== Find Algebraic Variables from ODE State Variables ===
==== Problem ====
When using lsode or other ODE solver to find the state variables x as a function of time t in a set of ODE's, there seems to be no simple way to return and tabulate algebraic variables y that are a function of the state variables x and time t, other than the derivatives.
Consider a modified version of the spring function defined above in the section on parametric functions:
<syntaxhighlight lang="Octave">
function sprime=dspring (s, t)
  k=1;
  tau=1;
  x = s(1);
  v = s(2);
  y=x+v*tau;
  sprime(1) = v;
  sprime(2) = -k * y;
endfunction
</syntaxhighlight>
What would be a way to way to use this function and lsode to tabulate not just the elements of x and sprime, but also (say) y as a function of time?
==== Solution ====
Modify the function sprime to return y as a second argument:
<syntaxhighlight lang="Octave">
function [sprime,y]=dspring (s, t)
  k=1;
  tau=1;
  x = s(1);
  v = s(2);
  y=x+v*tau;
  sprime(1) = v;
  sprime(2) = -k * y;
endfunction
</syntaxhighlight>
lsode will happily ignore the second return argument y, but after x has been tabulated, y can be tabulated using the second return argument of the spring function in a for loop, e.g.
<syntaxhighlight lang="Octave">
t = linspace (0, 10, 100);
x = lsode ('dspring',[1;0], t);
y=zeros(length(t),1);
for i=1:length(t)
  [xtmp,y(i,1)]=dspring(x(i,:),t(i));
end
</syntaxhighlight>
==== Discussion ====
This is a simple example with one algebraic variable. Others can be added by extending the second dimension of the y array.


[[Category:Tutorials]]
[[Category:Examples]]

Latest revision as of 23:58, 1 February 2024

An Octave cookbook. Each entry should go in a separate section and have the following subsection: problem, solution, discussion and maybe a see also.


Programs, Libraries, and Packages[edit]

Recipes for developers of Octave programs and libraries. The type of stuff an Octave programmer should be aware when writing code for others.

Find Octave configuration[edit]

Octave can be built with many configurations so programs may end up running in a machine without features they need. Developers should never expect an Octave installation to have all features. And programs should identify if the required features are available.

This is a list of possible tests to check for features:

## support for 64 bit indexing
sizemax () > intmax ("int32")

## built with support for java
usejava ("jvm")

## Image IO with support for tif files
any (cellfun (@(x) ismember ("tif", x), {imformats.ext}))
## Image IO with support for png files
any (cellfun (@(x) ismember ("png", x), {imformats.ext}))

Find if a package is installed[edit]

Problem[edit]

You have a program that uses different functions or behaves different depending on the availability of specific packages.

Solution[edit]

Use pkg ("list", pkg-name) like so:

if (! isempty (pkg ("list", "foo")))
  ## use functions from package foo, the preferred way
elseif (! isempty (pkg ("list", "bar")))
  ## use functions from package bar, not so optimal
else
  ## default case
endif

Discussion[edit]

It's not recommended to use this if the only purpose is to then fail in the absence of the package. In such case, simply try to load the package and Octave will already give a error message that is informative enough.

There is only purpose to check this, if there is something different to do if a package is missing. The same is true for catching an error from pkg load. If you only catch an error to then throw it again then you might as well not catch it in the first place.

## This contraption doesn't add anything.  If 'pkg load' fails, it
## will already give an error message to the user.
try
  pkg load foo;
catch
  error ("failed to load foo: %s", lasterr ());
end_try_catch

## Again, doesn't add anything.  The failure of 'pkg load' is enough
if (isempty (pkg ("list", "foo")))
  error ("program: package foo is not installed");
endif

Beware that an installed package is not always a guarantee that a function will be available. Some packages may disable functions at build time, or specific functions may have specific runtime requirements.

List all function in Octave[edit]

Use the following script (filename list_func.m)

## List of all builtin (C++) functions and m-file functions
funcs = vertcat (__builtins__ (), __list_functions__ ());

## Write list to file
fid = fopen ("all_funcs.tmp", "w");
if (fid == -1)
  error ("Unable to open temporary file all_funcs.tmp.  Aborting...\n");
endif
fprintf (fid, "%s\n", funcs{:});
fclose (fid);

And execute with

   run-octave -f list_func.m

Structures[edit]

Retrieve a field value from all entries in a struct array[edit]

Problem[edit]

You have a struct array with multiple fields, and you want to access the value from a specific field from all elements. For example, you want to return the age from all patients in the following case:

samples = struct ("patient", {"Bob", "Kevin", "Bob" , "Andrew"},
                  "age",     { 45  ,  52    ,  45   ,  23     },
                  "protein", {"H2B", "CDK2" , "CDK2", "Tip60" },
                  "tube"   , { 3   ,  5     ,  2    ,  18     }
                  );

Solution[edit]

Indexing the struct returns a comma separated list so use them to create a matrix.

[samples(:).age]

This however does not keep the original structure of the data, instead returning all values in a single column. To fix this, use reshape().

reshape ([samples(:).age], size (samples))

Discussion[edit]

Returning all values in a comma separated lists allows you to make anything out of them. If numbers are expected, create a matrix by enclosing them in square brackets. But if strings are to be expected, a cell array can also be easily generated with curly brackets

{samples(:).patient}

You are also not limited to return all elements, you may use logical indexing from other fields to get values from the others:

[samples([samples(:).age] > 34).tube]                 ## return tube numbers from all samples from patients older than 34
[samples(strcmp({samples(:).protein}, "CDK2")).tube]  ## return all tube numbers for protein CDK2

Array manipulation[edit]

Select a slice from an n-D array[edit]

Problem[edit]

For an array A with arbitrary number of dimensions, select, for example, the first column. This would be A(:, 1) if A was 2-D, A(:, 1, :) if A was 3-D, and so on.

Solution[edit]

One possibility is to use subsref with the input idx created dynamically with repelems to have the right number of dimensions. This can be written as a function:

function [B]= array_slice (A,k,d)
#return the k-th slice (row, column...) of A, with d specifying the dimension to slice on
  idx.type = "()";
  idx.subs = repelems ({':'}, [1;ndims(A)]);
  idx.subs(d) = k;
  B = subsref (A,idx);
endfunction

#test cases
%!shared A
%! A=rand(2, 3);
%!assert (array_slice (A,1,2), A(:, 1))
%! A=rand(2, 3, 4);
%!assert (array_slice (A,2,1), A(2, :, :))
%! A=rand(2, 3, 4, 5);
%!assert (array_slice (A,1,2), A(:, 1, :, :))
%! A=rand(2, 3, 4, 5, 6);
%!assert (array_slice (A,2,3), A(:, :, 2, :, :))

To remove the singleton dimension d from the result B, use

B = reshape(B, [size(B)([1:d-1 d+1:end])]);

Input/output[edit]

Display matched elements from different arrays[edit]

Problem[edit]

You have two, or more, arrays with paired elements and want to print out a string about them. For example:

keys   = {"human",  "mouse", "chicken"};
values = [ 64        72       70      ];

and you want to display:

Calculated human genome GC content is 64%
Calculated mouse genome GC content is 72%
Calculated chicken genome GC content is 70%

Solution[edit]

Make a two rows cell array, with each paired data in a column and supply a cs-list to printf

values = num2cell (values);
new    = {keys{:}; values{:}};
printf ("Calculated %s genome GC content is %i%%\n", new{:})

or in a single line:

printf ("Calculated %s genome GC content is %i%%\n", {keys{:}; num2cell(values){:}}{:})

Discussion[edit]

printf and family do not accept cell arrays as values. However, they keep repeating the template given as long as it has enough arguments to keep going. As such, the trick is on supplying a cs-list of elements which can be done by using a cell array and index it with {}.

Since values are stored in column-major order, paired values need to be on the same column. A new row of data can then be added later with

new(end+1,:) = {"Andrew", "Bob", "Kevin"};

Note that normal brackets are now being used for indexing.

Swap values[edit]

If you want to exchange the value between two variables without creating a dummy one, you can simply do:

[b,a] = deal (a,b);

Collect all output arguments of a function[edit]

If you have a function that returns several values, e.g.

function [a b c]= myfunc ()
  [a,b,c] = deal (1,2,3); 
endfunction

and you want to collect them all into a single cell (similarly to Python's zip() function) you can do:

outargs = nthargout (1:3, @myfunc)

Create a text table with fprintf[edit]

Imagine that you want to create a text table with fprintf with 2 columns of 15 characters width and both right justified. How to do this thing?

That's easy:

If the variable Text is a cell array of strings (of length < 15) with two columns and a certain number of rows, simply type for the k-th row of Text

fprintf('%15.15s | %15.15s\n', Text{k,1}, Text{k,2});

The syntax %<n>.<m>s allocates n places to write chars and display the m first characters of the string to display.

Example:

Text = {"Hello", "World"};
fprintf('%15.15s | %15.15s\n', Text{1,1}, Text{1,2})
         Hello |           World

Load comma separated values (*.csv) files[edit]

  1. Using textread gets you a one column cell array. The original size can be restored by using the reshape function.
    A = textread("file.csv", "%d", "delimiter", ",");
    B = textread("file.csv", "%s", "delimiter", ",");
    inds = isnan(A);
    B(! inds) = num2cell (A(! inds))
    
  2. Another option is to use the function csvread. However, this function can't handle non-numerical data.
  3. The probably best option is to use the function csv2cell from the IO package. This function can read mixed-type (numerical and text) .csv files, allows to specify other field separators than a comma and other text protection characters (default: " double quote) and can skip header lines. If you have the IO package installed and loaded, type help csv2cell at the Octave prompt for more info.

Load XML files[edit]

Reading XML in octave can be achieved using the java library Apache Xerces.

It seems that the Matlab's xmlread is just a thin wrapper around the Apache Xerces library. One should note however, that Java functions have the working directory set to the working directory when octave starts and the working directory is not modified by a cd in octave. Matlab has the same behavior, as Java does not provide a way to change the current working directory (http://bugs.java.com/bugdatabase/view_bug.do?bug_id=4045688). To avoid any issues, it is thus better to use the absolute path to the XML file.

You need the jar files xercesImpl.jar and xml-apis.jar from e.g. https://xerces.apache.org/mirrors.cgi#binary (check for the latest version). Use javaaddpath to include these files:

javaaddpath ("/path/to/xerces-2_11_0/xercesImpl.jar");
javaaddpath ("/path/to/xerces-2_11_0/xml-apis.jar");

Example:

filename = "sample.xml";

## These three lines are equivalent to xDoc = xmlread(filename) in Matlab
parser = javaObject("org.apache.xerces.parsers.DOMParser");
parser.parse(filename); 
xDoc = parser.getDocument();

elem = xDoc.getElementsByTagName("data").item(0); ## get first data element 
data = elem.getFirstChild.getTextContent();       ## get text from child
att  = elem.getAttribute("att");                  ## get attribute named att
File: sample.xml
<root>
  <data att="1">hello</data>
</root>

Using variable strings in commands[edit]

For example, to plot data using a string variable as a legend:

  1. Static string as legend (simplest):
    x = linspace (-1, 3, 100);
    y = sin (x);
    legend = "-1;My data;";
    plot (x, y, legend);
    
  2. Variable string as legend (moderate):
    x = linspace (-1, 3, 100);
    y = sin (x);
    dataName = "My data";
    plot (x, y, sprintf("-1;%s;", dataName));
    
  3. Variable string as legend using eval (not as neat):
    legend = "My data";
    plot_command = ["plot (x, y, '-1;", legend, ";')"];
    eval (plot_command);
    

These same tricks are useful for reading and writing data files with unique names, etc.

Combinatorics[edit]

Combinations with string characters[edit]

Problem[edit]

You want to get all combinations of different letters but nchoosek only accepts numeric input.

Solution[edit]

Convert your string to numbers and then back to characters.

string = "Hello";
n = 4;
char (nchoosek (uint8 (string), n))

Discussion[edit]

A string in Octave is just a character matrix and can easily be converted to numeric form back and forth. Each character has an associated number (the asci function of the miscellaneous package displays a nicely formatted conversion table).

Permutations with repetition[edit]

Problem[edit]

You want to generate all possible permutations of a vector with repetition.

Solution[edit]

Use ndgrid:

[x, y, z] = ndgrid ([1, 2, 3, 4, 5]);
[x(:), y(:), z(:)]

Discussion[edit]

It is possible to expand the code above and make it work for any length of permutations.

cart  = nthargout ([1:n], @ndgrid, vector);
combs = cell2mat (cellfun (@(c) c(:), cart, "UniformOutput", false));

Mathematics[edit]

Test if a number is an integer[edit]

The simplest method is probably fix:

fix (x) == x

Find if a number is even/odd[edit]

Problem[edit]

You have a number, or an array or matrix of them, and want to know if any of them is an odd or even number, i.e., their parity.

Solution[edit]

Check the remainder of a division by two. If the remainder is zero, the number is even.

mod (value, 2)  ## 1 if odd, zero if even

Since mod accepts a matrix, the following can be done:

any  (mod (values, 2))  ## true if at least one number in values is even
all  (mod (values, 2))  ## true if all numbers in values are odd

any (!logical (mod (values, 2)))  ## true if at least one number in values is even
all (!logical (mod (values, 2)))  ## true if all numbers in values are even

Discussion[edit]

Since we are checking for the remainder of a division, the first choice would be to use rem. However, in the case of negative numbers mod will still return a positive number making it easier for comparisons. Another alternative is to use bitand (X, 1) or bitget (X, 1) but those are a bit slower.

Note that this solution applies to integers only. Non-integers such as 0.5 or 4.201 are neither even nor odd. If the source of the numbers are unknown, such as user input, some sort of checking should be applied for NaN, Inf, or non-integer values.

Parametrized Functions[edit]

Problem[edit]

One sometimes needs to define a family of functions depending on a set of parameters, e.g., where denote a the variables on which the function operates and are the parameters used to chose one specific element of the family of functions.

For example, let's say we need to compute the time evolution of the elongation of a spring for different values of the spring constant

Solution[edit]

We could solve the problem with the following code to solve the spring equation for different values of the spring constant:

t = linspace (0, 10, 100);
function sprime = spring (s, t, k)
  x = s(1);
  v = s(2);
  sprime(1) = v;
  sprime(2) = -k * x;
endfunction
k = 1;
x1 = lsode (@(x, t) spring (x, t, k), [1;0], t)(:, 1);
k = 2;
x2 = lsode (@(x, t) spring (x, t, k), [1;0], t)(:, 2);
plot (t, x1, t, x2)
legend ('x1', 'x2')

Solparfun.png

Discussion[edit]

In the above example, the function "sprime" represents a family of functions of the variables parametrized by the parameter .

The anonymous function

@(x, t) sprime (x, t, k)

is a function of only where the parameter is "frozen" to the value it has at the moment in the current scope.

Distance between points[edit]

Problem[edit]

Given a set of points in space we want to calculate the distance between all of them. Each point is described by its components . Assume that the points are saved in a matrix P with m rows (one for each point) and n columns, one for each component.

Solution[edit]

One way of proceeding is to use the broadcast properties of operators in GNU Octave. The square distance between the points can be calculated with the code:

[m, n] = size (P);
Dsq     = zeros (m); 
for i = 1:n
  Dsq += (P(:,i) - P(:,i)').^2;
endfor

This matrix is symmetric with zero diagonal.

Similarly the vectors pointing from one point to the another is

R = zeros (m, m, n); 
for i = 1:n
  R(:,:,i) = P(:,i) - P(:,i)';
endfor

The relation between Dsq and R is

Dsq = sumsq (R, 3);

Discussion[edit]

The calculation can be implemented using functions like cellfun and avoid the loop over components of the points. However in most cases we will have more points than components and the improvement, if any, will be minimal.

Another observation is that the matrix Dsq is symmetric and we could store only the lower or upper triangular part. To use this optimization in a practical way check the help of the functions vech and unvech (this one is in the Forge package general). Two functions that haven't seen the light yet are sub2ind_tril and ind2sub_tril (currently private functions in the Forge package mechanics) that are useful to index the elements of a vector constructed with the function vech. Each page (the third index) of the multidimensional array R is an anti-symmetric matrix and we could also save some memory by keeping only one of the triangular submatrices.

Check the Geometry package for many more distance functions (points, lines, polygons, etc.).

Find Algebraic Variables from ODE State Variables[edit]

Problem[edit]

When using lsode or other ODE solver to find the state variables x as a function of time t in a set of ODE's, there seems to be no simple way to return and tabulate algebraic variables y that are a function of the state variables x and time t, other than the derivatives.

Consider a modified version of the spring function defined above in the section on parametric functions:

function sprime=dspring (s, t)
  k=1;
  tau=1;
  x = s(1);
  v = s(2);
  y=x+v*tau;
  sprime(1) = v;
  sprime(2) = -k * y;
endfunction

What would be a way to way to use this function and lsode to tabulate not just the elements of x and sprime, but also (say) y as a function of time?

Solution[edit]

Modify the function sprime to return y as a second argument:

function [sprime,y]=dspring (s, t)
  k=1;
  tau=1;
  x = s(1);
  v = s(2);
  y=x+v*tau;
  sprime(1) = v;
  sprime(2) = -k * y;
endfunction

lsode will happily ignore the second return argument y, but after x has been tabulated, y can be tabulated using the second return argument of the spring function in a for loop, e.g.

t = linspace (0, 10, 100);
x = lsode ('dspring',[1;0], t);
y=zeros(length(t),1);
for i=1:length(t)
  [xtmp,y(i,1)]=dspring(x(i,:),t(i));
end


Discussion[edit]

This is a simple example with one algebraic variable. Others can be added by extending the second dimension of the y array.