Cookbook: Difference between revisions

Jump to navigation Jump to search
2,280 bytes added ,  22 August 2012
m (→‎Parametrized Functions: no need to break lines when using the Math extension)
Line 98: Line 98:
=== 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.
==== 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| -->
{{Code|Calculate square distance between points|<syntaxhighlight lang="octave" style="font-size:13px">
[N, dim] = size (P);
Dsq    = zeros (N);
for i = 1:dim
  Dsq += (P(:,i) - P(:,i)').^2;
endfor
</syntaxhighlight>}}
This matrix is symmetric with zero diagonal.
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">
R    = zeros (N,N,dim);
for i = 1:dim
  R(:,:,i) = P(:,i) - P(:,i)';
endfor
</syntaxhighlight>}}
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>}}
==== 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.


== Plotting ==
== Plotting ==
== User input ==
== User input ==
657

edits

Navigation menu