Latest revision |
Your text |
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 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 acess 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"}, |
| samples = struct ("patient", {"Bob", "Kevin", "Bob" , "Andrew"}, | | "age", { 45 , 52 , 45 , 23 }, |
| "age", { 45 , 52 , 45 , 23 },
| | "protein", {"H2B", "CDK2" , "CDK2", "Tip60" }, |
| "protein", {"H2B", "CDK2" , "CDK2", "Tip60" },
| | "tube" , { 3 , 5 , 2 , 18 } |
| "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] |
| [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)) |
| 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 |
|
| |
|
| Returning all values in a comma separated lists allows you to make anything out of them.
| | {samples(:).patient} |
| 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([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 |
| [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 ===
| |
|
| |
| ==== Problem ====
| |
|
| |
| 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:
| |
|
| |
| Calculated human genome GC content is 64%
| |
| Calculated mouse genome GC content is 72%
| |
| Calculated chicken genome GC content is 70%
| |
|
| |
| ==== 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:
| |
|
| |
| <syntaxhighlight lang="Octave">
| |
| printf ("Calculated %s genome GC content is %i%%\n", {keys{:}; num2cell(values){:}}{:})
| |
| </syntaxhighlight>
| |
|
| |
| ==== Discussion ====
| |
|
| |
| {{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"> | | {{Code|Swap values without dummy variable|<syntaxhighlight lang="octave" style="font-size:13px"> |
| [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"> | | {{Code||<syntaxhighlight lang="octave" style="font-size:13px"> |
| 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"> | | {{Code|Collect multiple output arguments|<syntaxhighlight lang="octave" style="font-size:13px"> |
| 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 odd. |
|
| |
|
| 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 |
|
| |
|
| <syntaxhighlight lang="Octave">
| | Since {{Codeline|mod()}} acceps a matrix, the following can be done: |
| mod (value, 2) ## 1 if odd, zero if even | |
| </syntaxhighlight>
| |
|
| |
|
| Since {{manual|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 |
| <syntaxhighlight lang="Octave">
| | |
| any (mod (values, 2)) ## true if at least one number in values is even | | any (!logical (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 | | all (!logical (mod (values, 2))) ## true if all numbers in values are even |
| | |
| 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. |
|
| |
|
| Since we are checking for the remainder of a division, the first choice would be to use {{manual|rem}}.
| | 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. |
| 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.
| | ==== See also ==== |
| Non-integers such as <code>0.5</code> or <code>4.201</code> are neither even nor odd.
| | Find if a number is an integer. |
| 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 459: |
Line 87: |
|
| |
|
| ==== Solution ==== | | ==== Solution ==== |
| | We could solve the problem with the following code: |
|
| |
|
| We could solve the problem with the following code to solve the spring equation for different values of the spring constant:
| | {{Code|Solve spring equation for different values of the spring constant|<syntaxhighlight lang="octave" style="font-size:13px"> |
| | |
| <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 476: |
Line 103: |
| 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 484: |
Line 111: |
| 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 anonymous function] | | The [http://www.gnu.org/software/octave/doc/interpreter/Anonymous-Functions.html#Anonymous-Functions anonympus function] |
| | | <pre> |
| <syntaxhighlight lang="Octave"> | | @(x, t) sprime (x, t, k) |
| @(x, t) sprime (x, t, k) | | </pre> |
| </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. | | 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 |
| One way of proceeding is to use the broadcast properties of operators in GNU Octave. | | <!-- {{SyntaxHighlight| --> |
| The square distance between the points can be calculated with the code: | | {{Code|Calculate square distance between points|<syntaxhighlight lang="octave" style="font-size:13px"> |
| | | [N, dim] = size (P); |
| <syntaxhighlight lang="Octave"> | | Dsq = zeros (N); |
| [m, n] = size (P); | | for i = 1:dim |
| 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| --> |
| <syntaxhighlight lang="Octave"> | | {{Code|Calculate radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px"> |
| R = zeros (m, m, n); | | R = zeros (N,N,dim); |
| for i = 1:n | | for i = 1:dim |
| 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| --> |
| | {{Code|Calculate square distance from radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px"> |
| | Dsq = sumsq (R,3); |
| | </syntaxhighlight>}} |
|
| |
|
| <syntaxhighlight lang="Octave">
| |
| 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. |
|
| |
|
| The calculation can be implemented using functions like {{manual|cellfun}} and avoid the loop over components of the points.
| | 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. |
| 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 === | | == Plotting == |
| | | == User input == |
| ==== 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:Examples]]
| |