Attached is a function to compute the thing in subject line (don't think
I can spell it correctly twice). The c.h.f. itself is pretty simple,
and I attached tests from A&S plus a test or two of special cases.
This thing will create tables of the c.h.f. for appropriately sized
arrays and in order to do that it requires a second function called
"tablify", also attached. This function, and I hate the name so
suggestions are welcome, is an extension of common_size and could be
replaced by bsxfun if bsxfun took multiple array arguments and was more
flexible about output arguments.
If you feel this is useful, it could go in the -forge special functions
package.
Bob
## Copyright (C) 2012 Robert T. Short
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{y},@var{p}]} = confluent
(@var{a},@var{b},@var{z},@var{tol}=@var{eps})
##
## Compute the confluent hypergeometric function
## M(@var{a},@var{b},@var{z})
## using the series representation.
## M is also called 1_F_1 or Phi in much of the literature.
##
## The arguments @var{a}, @var{b}, and @var{z} are described
## in, for example, Abramowitz and Stegun. The first two,
## @var{a} and @var{b}, are the parameters and @var{z} is the
## argument.
##
## The confluent hypergeometric function is computed
## using the series expression for M, truncated when
## the coefficient is less than @var{tol}.
##
## If the second return argument @var{p} is present
## the coefficients of the truncated series are returned
## in @var{p}. This is useful when @var{a} is a negative
## integer and M is a polynomial.
##
## If the arguments are vectors or matrices with appropriate
## dimensions a table will be made of the results and p will
## be a cell array of polynomials.
##
## See also: @seealso {tablify}
##
## Reference: Abramowitz and Stegun
##
## @end deftypefn
## Author: Robert T. Short (rtsh...@ieee.org)
## Description: Confluent hypergeometric function.
function [y,p] = confluent(a,b,z,tol=eps)
if ( (nargin < 3) || (nargin > 4) )
usage("y = confluent(a,b,z,tol=eps)");
end
% Make a table and call the subfunction chf
% for each element in the table.
[a,b,z] = tablify(a,b,z);
y = zeros(size(a));
if (nargout>1)
p = cell(size(a));
end
for idx = 1:rows(a)
for jdx = 1:columns(a)
if (nargout>1)
[y(idx,jdx),p{idx,jdx}] = chf(a(idx,jdx),b(idx,jdx),z(idx,jdx),tol);
else
[y(idx,jdx)] = chf(a(idx,jdx),b(idx,jdx),z(idx,jdx),tol);
end
end
end
end
% Subfunction that actually computes the chf.
function [y,p] = chf(a,b,z,tol)
t = 1;
y = 1;
k = 0;
x = 1;
c = 1;
p = [];
do
if (nargout > 1)
p = [c, p];
end
k++;
c = c * (a+k-1)/(b+k-1)/k;
x = x*z;
t = c*x;
y = y + c*x;
until ( (abs(t)<tol) || (b+k-1==0) )
end
% Tests from Abramowitz and Stegun.
%
%!test
%! a = [-1.0;-0.5;0.0;+0.5;+1.0];
%! b = [ 0.1, 0.5, 1.0 ];
%! x = 0.1;
%! ct = [ [ 0.0000000 8.0000000e-01 9.0000000e-01 ];
%! [ 4.8836025e-01 8.9829940e-01 9.4936442e-01 ];
%! [ 1.0000000e+00 1.0000000e+00 1.0000000e+00 ];
%! [ 1.5354828e+00 1.1051709e+00 1.0519282e+00 ];
%! [ 2.0953812e+00 1.2138822e+00 1.1051709e+00 ]];
%! assert(confluent(a,b,x),ct,-5e-8);
%!
%!test
%! a = [-1.0;-0.5;0.0;+0.5;+1.0];
%! b = [ 0.1, 0.5, 1.0 ];
%! x = 2.0;
%! ct = [ [ -1.9000000e+01 -3.0000000e+00 -1.0000000e+00 ];
%! [ -1.6997468e+01 -2.0687595e+00 -3.6900042e-01 ];
%! [ 1.0000000e+00 1.0000000e+00 1.0000000e+00 ];
%! [ 4.5608543e+01 7.3890561e+00 3.4415239e+00 ];
%! [ 1.3143241e+02 1.8678878e+01 7.3890561e+00 ]];
%! assert(confluent(a,b,x),ct,-5e-8);
%!
%!test
%! a = [-1.0;-0.5;0.0;+0.5;+1.0];
%! b = [ 0.1, 0.5, 1.0 ];
%! x = 4.0;
%! ct = [ [ -3.9000000e+01 -7.0000000e+00 -3.0000000e+00 ];
%! [ -8.8670480e+01 -1.1212361e+01 -3.5187312e+00 ];
%! [ 1.0000000e+00 1.0000000e+00 1.0000000e+00 ];
%! [ 4.8014767e+02 5.4598150e+01 1.6843984e+01 ];
%! [ 1.8088849e+03 1.9364005e+02 5.4598150e+01 ]];
%! assert(confluent(a,b,x),ct,-5e-8);
%!
%!test
%! a = [-1.0;-0.5;0.0;+0.5;+1.0];
%! b = [ 0.1, 0.5, 1.0 ];
%! x = 6.0;
%! ct = [ [ -5.9000000e+01 -1.1000000e+01 -5.0000000e+00 ];
%! [ -4.6224363e+02 -4.9331877e+01 -1.3733318e+01 ];
%! [ 1.0000000e+00 1.0000000e+00 1.0000000e+00 ];
%! [ 4.2706845e+03 4.0342879e+02 9.8033340e+01 ];
%! [ 1.9250691e+04 1.7515977e+03 4.0342879e+02 ]];
%! assert(confluent(a,b,x),ct,-5e-8);
%!
%!test
%! a = [-1.0;-0.5;0.0;+0.5;+1.0];
%! b = [ 0.1, 0.5, 1.0 ];
%! x = 8.0;
%! ct = [ [ -7.9000000e+01 -1.5000000e+01 -7.0000000e+00 ];
%! [ -2.6742961e+03 -2.4320200e+02 -5.6658271e+01 ];
%! [ 1.0000000e+00 1.0000000e+00 1.0000000e+00 ];
%! [ 3.5774528e+04 2.9809580e+03 6.1706403e+02 ];
%! [ 1.8427980e+05 1.4944361e+04 2.9809580e+03 ]];
%! assert(confluent(a,b,x),ct,-5e-8);
%!
%!test
%! a = [-1.0;-0.5;0.0;+0.5;+1.0];
%! b = [ 0.1, 0.5, 1.0 ];
%! x = 10.0;
%! ct = [ [ -9.9000000e+01 -1.9000000e+01 -9.0000000e+00 ];
%! [ -1.6608619e+04 -1.3381435e+03 -2.6750359e+02 ];
%! [ 1.0000000e+00 1.0000000e+00 1.0000000e+00 ];
%! [ 2.9071300e+05 2.2026466e+04 4.0427554e+03 ];
%! [ 1.6645066e+06 1.2345819e+05 2.2026466e+04 ]];
%! assert(confluent(a,b,x),ct,-5e-8);
%
% Laguerre polynomials.
%
%!test
%! n=3;alpha=0;
%! [~,p] = confluent(-n,alpha+1,1);
%! p = p{1,1};
%! p = gamma(n+alpha+1)/gamma(alpha+1) * p;
%! assert(p, [-1,9,-18,6], eps);
## Copyright (C) 2012 Robert T. Short
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{y1}, @dots{}] =} tablify (@var{x1},
@dots{})
## Create a table out of the input arguments, if possible. The table is created
## by extending row and column vectors to like dimensions. If the dimensions of
## input vectors are not commensurate an error will occur. Dimensions are
## commensurate if they have the same number of rows or a single row and the
## same number of columns or a single column.
##
## For example:
##
## @example
## @group
## [a, b] = tablify ([1 2; 3 4], 5)
## @result{} a = [ 1, 2; 3, 4 ]
## @result{} b = [ 5, 5; 5, 5 ]
## @end group
## @end example
## @example
## @group
## [a, b, c] = tablify (1, [1 2 3 4], [5;6;7])
## @result{} b = [ 1 1 1 1; 1 1 1 1; 1 1 1 1]
## @result{} b = [ 1 2 3 4; 1 2 3 4; 1 2 3 4]
## @result{} c = [ 5 5 5 5; 6 6 6 6; 7 7 7 7 ]
## @end group
## @end example
##
## @noindent
## The following example attempts to expand vectors that do not
## have commensurate dimensions and will produce an error.
##
## @example
## @group
## tablify([1 2],[3 4 5])
## @end group
## @end example
##
## @noindent
## Note that use of array operations and broadcasting is more efficient
## for many situations.
##
## @seealso {common_size}
##
## @end deftypefn
## Author: Robert T. Short
## Created: 3/6/2012
function [varargout] = tablify (varargin)
if (nargin < 2)
varargout = varargin;
return;
endif
empty = cellfun (@isempty, varargin);
nrows = cellfun (@rows, varargin(~empty));
ridx = nrows>1;
if (any(ridx))
rdim = nrows(ridx)(1);
else
rdim = 1;
endif
ncols = cellfun (@columns, varargin(~empty));
cidx = ncols>1;
if (any(cidx))
cdim = ncols(cidx)(1);
else
cdim = 1;
endif
if ( any(nrows(ridx) != rdim ) || any(ncols(cidx) != cdim ))
error('tablify: incommensurate sizes');
endif
varargout = varargin;
varargout(~ridx) = cellindexmat(varargout(~ridx),ones(rdim,1),':');
varargout(~cidx) = cellindexmat(varargout(~cidx),':',ones(1,cdim));
endfunction
%!error tablify([1,2],[3,4,5])
%!error tablify([1;2],[3;4;5])
%!test
%! a = 1; b = 1; c = 1;
%! assert(tablify(a),a);
%! [x,y,z]=tablify(a,b,c);
%! assert(x,a);
%! assert(y,b);
%! assert(z,c);
%!test
%! a = 1; b = [1 2 3];
%! [x,y]=tablify(a,b);
%! assert(x,[1 1 1]);
%! assert(y,[1 2 3]);
%!test
%! a = 1; b = [1;2;3];
%! [x,y]=tablify(a,b);
%! assert(x,[1;1;1]);
%! assert(y,[1;2;3]);
%!test
%! a = [1 2]; b = [1;2;3]; c=1;
%! [x,y,z]=tablify(a,b,c);
%! assert(x,[1 2; 1 2; 1 2]);
%! assert(y,[1 1; 2 2; 3 3]);
%! assert(z,[1 1; 1 1; 1 1]);
%!test
%! a = [1 2]; b = [1;2;3]; c=[2 3];
%! [x,y,z]=tablify(a,b,c);
%! assert(x,[1 2; 1 2; 1 2]);
%! assert(y,[1 1; 2 2; 3 3]);
%! assert(z,[2 3; 2 3; 2 3]);
%!test
%! a = [1 2]; b = [1;2;3]; c=[];
%! [x,y,z]=tablify(a,b,c);
%! assert(x,[1 2; 1 2; 1 2]);
%! assert(y,[1 1; 2 2; 3 3]);
%! assert(z,[]);
------------------------------------------------------------------------------
This SF email is sponsosred by:
Try Windows Azure free for 90 days Click Here
http://p.sf.net/sfu/sfd2d-msazure
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
