2
edits
(32 intermediate revisions by 13 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 == | |||
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 === | |||
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: | |||
<syntaxhighlight lang="Octave"> | |||
## 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})) | |||
</syntaxhighlight> | |||
=== Find if a package is installed === | |||
==== Problem ==== | |||
You have a program that uses different functions or behaves different | |||
depending on the availability of specific packages. | |||
==== Solution ==== | |||
Use {{codeline|pkg ("list", pkg-name)}} like so: | |||
<syntaxhighlight lang="Octave"> | |||
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 | |||
</syntaxhighlight> | |||
==== Discussion ==== | |||
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 | |||
{{codeline|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. | |||
<syntaxhighlight lang="Octave"> | |||
## 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 | |||
</syntaxhighlight> | |||
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 === | |||
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 | 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: | ||
<syntaxhighlight lang="Octave"> | |||
samples = struct ("patient", {"Bob", "Kevin", "Bob" , "Andrew"}, | |||
"age", { 45 , 52 , 45 , 23 }, | |||
"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. | ||
<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()}}. | ||
<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 | |||
<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: | ||
<syntaxhighlight lang="Octave"> | |||
[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 == | |||
=== Select a slice from an n-D array === | |||
==== 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. | |||
==== Solution ==== | |||
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: | |||
<syntaxhighlight lang="Octave"> | |||
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, :, :)) | |||
</syntaxhighlight> | |||
To remove the singleton dimension {{Codeline|d}} from the result {{Codeline|B}}, use | |||
<syntaxhighlight lang="Octave"> | |||
B = reshape(B, [size(B)([1:d-1 d+1:end])]); | |||
</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: | ||
<syntaxhighlight lang="Octave"> | |||
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 mouse genome GC content is 72% | |||
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 {{manual|printf}} | |||
<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: | ||
<syntaxhighlight lang="Octave"> | |||
printf ("Calculated %s genome GC content is %i%%\n", {keys{:}; num2cell(values){:}}{:}) | |||
</syntaxhighlight> | |||
==== Discussion ==== | ==== Discussion ==== | ||
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 | {{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: | ||
<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. | ||
<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: | ||
<syntaxhighlight lang="Octave"> | |||
outargs = nthargout (1:3, @myfunc) | outargs = nthargout (1:3, @myfunc) | ||
</syntaxhighlight> | |||
=== Create a text table with fprintf=== | |||
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: | |||
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 | |||
<syntaxhighlight lang="Octave"> | |||
fprintf('%15.15s | %15.15s\n', Text{k,1}, Text{k,2}); | |||
</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: | |||
<syntaxhighlight lang="Octave"> | |||
Text = {"Hello", "World"}; | |||
fprintf('%15.15s | %15.15s\n', Text{1,1}, Text{1,2}) | |||
</syntaxhighlight> | |||
Hello | World | |||
===Load comma separated values (*.csv) files=== | |||
# 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", ","); | |||
B = textread("file.csv", "%s", "delimiter", ","); | |||
inds = isnan(A); | |||
B(! inds) = num2cell (A(! inds)) | |||
</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. | |||
===Load XML files=== | |||
Reading XML in octave can be achieved using the java library [https://xerces.apache.org/ Apache Xerces]. | |||
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. | |||
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: | |||
<syntaxhighlight lang="Octave"> | |||
javaaddpath ("/path/to/xerces-2_11_0/xercesImpl.jar"); | |||
javaaddpath ("/path/to/xerces-2_11_0/xml-apis.jar"); | |||
</syntaxhighlight> | |||
Example: | |||
<syntaxhighlight lang="Octave"> | |||
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 | |||
</syntaxhighlight> | |||
{{File|sample.xml| | |||
<syntaxhighlight lang="xml"> | |||
<root> | |||
<data att="1">hello</data> | |||
</root> | |||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
===Using variable strings in commands=== | |||
For example, to plot data using a string variable as a legend: | |||
# Static string as legend (simplest): <syntaxhighlight lang="Octave"> | |||
x = linspace (-1, 3, 100); | |||
y = sin (x); | |||
legend = "-1;My data;"; | |||
plot (x, y, legend); | |||
</syntaxhighlight> | |||
# Variable string as legend (moderate): <syntaxhighlight lang="Octave"> | |||
x = linspace (-1, 3, 100); | |||
y = sin (x); | |||
dataName = "My data"; | |||
plot (x, y, sprintf("-1;%s;", dataName)); | |||
</syntaxhighlight> | |||
# Variable string as legend using {{manual|eval}} (not as neat): <syntaxhighlight lang="Octave"> | |||
legend = "My data"; | |||
plot_command = ["plot (x, y, '-1;", legend, ";')"]; | |||
eval (plot_command); | |||
</syntaxhighlight> | |||
These same tricks are useful for reading and writing data files with unique names, etc. | |||
== Combinatorics == | |||
=== Combinations with string characters === | |||
==== Problem ==== | |||
You want to get all combinations of different letters but {{manual|nchoosek}} only accepts numeric input. | |||
==== Solution ==== | |||
Convert your string to numbers and then back to characters. | |||
<syntaxhighlight lang="Octave"> | |||
string = "Hello"; | |||
n = 4; | |||
char (nchoosek (uint8 (string), n)) | |||
</syntaxhighlight> | |||
==== 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). | |||
=== Permutations with repetition === | |||
==== Problem ==== | |||
You want to generate all possible permutations of a vector with repetition. | |||
==== Solution ==== | |||
Use {{manual|ndgrid}}: | |||
<syntaxhighlight lang="Octave"> | |||
[x, y, z] = ndgrid ([1, 2, 3, 4, 5]); | |||
[x(:), y(:), z(:)] | |||
</syntaxhighlight> | |||
==== Discussion ==== | |||
It is possible to expand the code above and make it work for any length of permutations. | |||
<syntaxhighlight lang="Octave"> | |||
cart = nthargout ([1:n], @ndgrid, vector); | |||
combs = cell2mat (cellfun (@(c) c(:), cart, "UniformOutput", false)); | |||
</syntaxhighlight> | |||
== Mathematics == | == Mathematics == | ||
=== 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 even. | |||
<syntaxhighlight lang="Octave"> | |||
mod (value, 2) ## 1 if odd, zero if even | |||
</syntaxhighlight> | |||
Since {{manual|mod}} accepts a matrix, the following can be done: | |||
<syntaxhighlight lang="Octave"> | |||
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 | |||
</syntaxhighlight> | |||
==== Discussion ==== | ==== Discussion ==== | ||
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. | |||
Note that this solution applies to integers only. | |||
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 116: | Line 459: | ||
==== Solution ==== | ==== Solution ==== | ||
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 132: | 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 140: | Line 484: | ||
In the above example, the function "sprime" represents a family of functions of the variables <math>x, t</math> parametrized by the parameter <math>k</math>. | In the above example, the function "sprime" represents a family of functions of the variables <math>x, t</math> parametrized by the parameter <math>k</math>. | ||
The [http://www.gnu.org/software/octave/doc/interpreter/Anonymous-Functions.html#Anonymous-Functions | The [http://www.gnu.org/software/octave/doc/interpreter/Anonymous-Functions.html#Anonymous-Functions anonymous function] | ||
is a function of only <math>x, t</math> where the parameter <math>k</math> is | <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>. | |||
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 | |||
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: | |||
[ | |||
Dsq = zeros ( | <syntaxhighlight lang="Octave"> | ||
for i = 1: | [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 lang="Octave"> | |||
R | R = zeros (m, m, n); | ||
for i = 1: | for i = 1:n | ||
R(:,:,i) = P(:,i) - P(:,i)'; | R(:,:,i) = P(:,i) - P(:,i)'; | ||
endfor | endfor | ||
</syntaxhighlight> | </syntaxhighlight> | ||
The relation between < | The relation between <code>Dsq</code> and <code>R</code> is | ||
<syntaxhighlight lang="Octave"> | |||
Dsq = sumsq (R, 3); | |||
</syntaxhighlight> | |||
==== Discussion ==== | ==== Discussion ==== | ||
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 < | 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) | |||
</syntaxhighlight> | [xtmp,y(i,1)]=dspring(x(i,:),t(i)); | ||
end | |||
</syntaxhighlight> | |||
==== Discussion ==== | |||
This | This is a simple example with one algebraic variable. Others can be added by extending the second dimension of the y array. | ||
[[Category:Examples]] |
edits