Mercurial > octave

changeset 24909:4a341330ee15 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added sine integral and cosine integral functions. -- added scripts/specfun/cosint.m added scripts/specfun/sinint.m changed NEWS changed doc/interpreter/arith.txi changed scripts/specfun/expint.m changed scripts/specfun/module.mk
author Michele Ginesi <michele.ginesi@gmail.com>
date Thu, 07 Sep 2017 16:28:35 +0200
parents 6c082a43abd8
children b98755ef7572
files NEWS doc/interpreter/arith.txi scripts/specfun/cosint.m scripts/specfun/expint.m scripts/specfun/module.mk scripts/specfun/sinint.m
diffstat 6 files changed, 394 insertions(+), 2 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS	Thu Sep 07 16:24:15 2017 +0200
+++ b/NEWS	Thu Sep 07 16:28:35 2017 +0200
@@ -252,6 +252,7 @@
       camva
       camzoom
       corrcoef
+      cosint
       erase
       gammaincinv
       getframe
@@ -270,6 +271,7 @@
       repelem
       rgb2gray
       rticks
+      sinint
       tfqmr
       thetaticks
       vecnorm
--- a/doc/interpreter/arith.txi	Thu Sep 07 16:24:15 2017 +0200
+++ b/doc/interpreter/arith.txi	Thu Sep 07 16:28:35 2017 +0200
@@ -295,6 +295,8 @@
 
 @DOCSTRING(commutation_matrix)
 
+@DOCSTRING(cosint)
+
 @DOCSTRING(duplication_matrix)
 
 @DOCSTRING(dawson)
@@ -330,6 +332,8 @@
 
 @DOCSTRING(psi)
 
+@DOCSTRING(sinint)
+
 @node Rational Approximations
 @section Rational Approximations
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/specfun/cosint.m	Thu Sep 07 16:28:35 2017 +0200
@@ -0,0 +1,200 @@
+## Copyright (C) 2017 Michele Ginesi
+##
+## 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/>.
+
+## Author: Michele Ginesi <michele.ginesi@gmail.com>
+
+## -*- texinfo -*-
+## @deftypefn {} {} cosint (@var{x})
+## Compute the cosine integral function:
+## @tex
+## $$
+## {\rm Ci} (x) = - \int_x^\infty {{\cos (t)} \over t} dt
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##             +oo
+##            /
+## Ci (x) = - | (cos (t)) / t dt
+##            /
+##           x
+##
+## @end group
+## @end example
+## @end ifnottex
+## An equivalent definition is
+## @tex
+## $$
+## {\rm Ci} (x) = \gamma + \log (x) + \int_0^x {{\cos (t) - 1} \over t} dt
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##
+##                              x
+##                             /
+##                             |   cos (t) - 1
+## Ci (x) = gamma + log (x) +  | ---------------  dt
+##                             |         t
+##                             /
+##                            0
+## @end group
+## @end example
+## @end ifnottex
+## Reference:
+##
+## @nospell{M. Abramowitz and I.A. Stegun},
+## @cite{Handbook of Mathematical Functions}
+## 1964.
+##
+## @seealso{expint, cos, sinint}
+##
+## @end deftypefn
+
+function [y] = cosint (x)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  sz = size (x);
+  x = x(:);
+  y = zeros (size (x));
+
+  i_miss = true (length (x), 1);
+
+  # Trivial values
+  y(x == 0) = - Inf;
+  y(x == Inf) = 0;
+  y(x == - Inf) = 1i * pi;
+
+  i_miss = ((i_miss) & (x != 0) & (x != Inf) & (x != - Inf));
+
+  # For values large in module and not in (-oo,0),we use the relation
+  # with expint
+
+  flag_large = (abs (x) > 2 & ((abs (imag (x)) > 1e-15) | real (x) > 0));
+  xx = x(flag_large);
+
+  # Abramowitz, relation 5.2.20
+  ii_sw = (real (xx) <= 0 & imag (xx) <= 0);
+  xx(ii_sw) = conj (xx(ii_sw));
+  ii_nw = (real (xx) < 0);
+  xx(ii_nw) *= -1;
+  yy = -0.5 * (expint (1i * xx) + expint (-1i * xx));
+  yy(ii_nw) += 1i * pi;
+  yy(ii_sw) = conj (yy(ii_sw));
+  y(i_miss & flag_large) = yy;
+
+  # For values small in module we use the series expansion
+
+  i_miss = ((i_miss) & (!flag_large));
+  xx = x(i_miss);
+  ssum = - xx .^ 2 / 4; # First term of the series expansion
+  gma = 0.57721566490153286060651209008; # Euler gamma constant
+  yy = gma + log (xx) + ssum;
+  flag_sum = true (nnz (i_miss), 1);
+  it = 1;
+  maxit = 300;
+  while ((any (flag_sum)) && (it < maxit));
+    ssum .*= - xx .^ 2 * (2 * it) / ((2 * it + 2) ^ 2 * (2 * it + 1));
+    yy(flag_sum) += ssum (flag_sum);
+    flag_sum = (abs (ssum) >= eps);
+    it++;
+  endwhile
+
+  y(i_miss) = yy;
+
+  y = reshape (y, sz);
+
+endfunction
+
+%!assert (cosint (1.1), 0.38487337742465081550, 2 * eps);
+
+
+%!test
+%! x = [2, 3, pi; exp(1), 5, 6];
+%! A = cosint (x);
+%! B = [0.422980828774864996, 0.119629786008000328, 0.0736679120464254860; ...
+%!      0.213958001340379779, -0.190029749656643879, -0.0680572438932471262];
+%! assert (A, B, -5e-15);
+
+%!assert (cosint (0), - Inf)
+%!assert (cosint (inf), 0)
+%!assert (cosint (-inf), 1i * pi)
+
+%%tests against maple
+%!assert (cosint (1), 0.337403922900968135, -2*eps)
+%!assert (cosint (-1), 0.337403922900968135 + 3.14159265358979324*I, -2*eps)
+%!assert (cosint (pi), 0.0736679120464254860, -2*eps)
+%!assert (cosint (-pi), 0.0736679120464254860 + 3.14159265358979324*I, -2*eps)
+%!assert (cosint (300), -0.00333219991859211178, -2*eps)
+%!assert (cosint (1e4), -0.0000305519167244852127, -2*eps)
+%!assert (cosint (20i), 1.28078263320282944e7 + 1.57079632679489662*I, -2*eps)
+
+%!test
+%! x = (0:4).';
+%! y_ex = [-Inf
+%!         0.337403922900968135
+%!         0.422980828774864996
+%!         0.119629786008000328
+%!         -0.140981697886930412];
+%! assert (cosint(x), y_ex, -3e-15);
+
+%!test
+%! x = -(0:4) ';
+%! y_ex = [-Inf
+%!         0.337403922900968135 + 3.14159265358979324*I
+%!         0.422980828774864996 + 3.14159265358979324*I
+%!         0.119629786008000328 + 3.14159265358979324*I
+%!         -0.140981697886930412 + 3.14159265358979324*I];
+%! assert (cosint (x), y_ex, -4*eps);
+
+%!test
+%! x = 1i * (0:4).';
+%! y_ex = [-Inf
+%!         0.837866940980208241 + 1.57079632679489662*I
+%!         2.45266692264691452 + 1.57079632679489662*I
+%!         4.96039209476560976 + 1.57079632679489662*I
+%!         9.81354755882318556 + 1.57079632679489662*I];
+%! assert (cosint (x), y_ex, -4*eps);
+
+%!test
+%! x = - 1i * (0:4).';
+%! y_ex = [- Inf
+%!         0.837866940980208241 - 1.57079632679489662*I
+%!         2.45266692264691452 - 1.57079632679489662*I
+%!         4.96039209476560976 - 1.57079632679489662*I
+%!         9.81354755882318556 - 1.57079632679489662*I];
+%! assert (cosint (x), y_ex, -4*eps);
+
+%!test
+%! x = [1+2i; -2+5i; 2-5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
+%! A = [ 2.03029639329172164 - 0.151907155175856884*I
+%!      1.61538963829107749 + 19.7257540553382650*I
+%!      1.61538963829107749 + 16.5841614017484717*I
+%!      -0.00514882514261049214
+%!      1246.11448604245441 + 1.57079632679489662*I
+%!      -8.63307471207423322 + 3.13159298695312800*I
+%!      0.0698222284673061493 - 3.11847446254772946*I ];
+%! B = cosint (x);
+%! assert (A, B, -6e-16)
--- a/scripts/specfun/expint.m	Thu Sep 07 16:24:15 2017 +0200
+++ b/scripts/specfun/expint.m	Thu Sep 07 16:28:35 2017 +0200
@@ -86,7 +86,7 @@
 ## @cite{Asymptotic expansions of integrals}
 ## 1986.
 ##
-## @seealso{exp}
+## @seealso{cosint, exp, sinint}
 ##
 ## @end deftypefn
 
--- a/scripts/specfun/module.mk	Thu Sep 07 16:24:15 2017 +0200
+++ b/scripts/specfun/module.mk	Thu Sep 07 16:28:35 2017 +0200
@@ -5,6 +5,7 @@
   %reldir%/betainc.m \
   %reldir%/betaincinv.m \
   %reldir%/betaln.m \
+  %reldir%/cosint.m \
   %reldir%/ellipke.m \
   %reldir%/expint.m \
   %reldir%/factor.m \
@@ -21,7 +22,8 @@
   %reldir%/primes.m \
   %reldir%/reallog.m \
   %reldir%/realpow.m \
-  %reldir%/realsqrt.m
+  %reldir%/realsqrt.m \
+  %reldir%/sinint.m
 
 %canon_reldir%dir = $(fcnfiledir)/specfun
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/specfun/sinint.m	Thu Sep 07 16:28:35 2017 +0200
@@ -0,0 +1,184 @@
+## Copyright (C) 2017 Michele Ginesi
+##
+## 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/>.
+
+## Author: Michele Ginesi <michele.ginesi@gmail.com>
+
+## -*- texinfo -*-
+## @deftypefn {} {} sinint (@var{x})
+## Compute the sine integral function:
+## @tex
+## $$
+## {\rm Si} (x) = \int_0^x {\sin (t) \over t} dt
+## $$
+## @end tex
+## @ifnottex
+##
+## @example
+## @group
+##            x
+##           /
+## Si (x) =  | sin (t) / t dt
+##           /
+##          0
+## @end group
+## @end example
+## @end ifnottex
+##
+## Reference:
+##
+## @nospell{M. Abramowitz and I.A. Stegun},
+## @cite{Handbook of Mathematical Functions}
+## 1964.
+##
+## @seealso{cosint, expint, sin}
+##
+## @end deftypefn
+
+function [y] = sinint (x)
+
+  if (nargin != 1)
+    print_usage ();
+  endif
+
+  sz = size (x);
+  x = x(:);
+  y = zeros (size (x));
+
+  i_miss = true (length (x), 1);
+
+  # Trivial values
+  y(x == 0) = 0;
+  y(x == Inf) = pi / 2;
+  y(x == - Inf) = - pi / 2;
+
+  i_miss = ((i_miss) & (x != 0) & (x != Inf) & (x != - Inf));
+
+  # For values large in module we use the relation with expint
+
+  flag_large = abs (x) > 2;
+  xx = x(flag_large & i_miss);
+  ii_neg = (real (xx) < 0);
+  xx(ii_neg) *= -1;
+  ii_conj = ((real (xx) == 0) & (imag (xx) < 0));
+  xx(ii_conj) = conj (xx(ii_conj));
+  yy = -0.5i * (expint (1i * xx) - expint (-1i * xx)) + pi / 2;
+  yy(ii_neg) *= -1;
+  yy(ii_conj) = conj (yy(ii_conj));
+  y(i_miss & flag_large) = yy;
+
+  # For values small in module we use the series expansion
+
+  i_miss = ((i_miss) & (!flag_large));
+  xx = x(i_miss);
+  ssum = xx; # First term of the series expansion
+  yy = ssum;
+  flag_sum = true (nnz (i_miss), 1);
+  it = 0;
+  maxit = 300;
+  while ((any (flag_sum)) && (it < maxit));
+    ssum .*= - xx .^ 2 * (2 * it + 1) / ((2 * it + 3) ^ 2 * (2 * it + 2));
+    yy(flag_sum) += ssum (flag_sum);
+    flag_sum = (abs (ssum) >= eps);
+    it++;
+  endwhile
+
+  y(i_miss) = yy;
+
+  y = reshape (y, sz);
+
+endfunction
+
+%!test
+%! x = 1.1;
+%! %y = sym(11)/10;
+%! A = sinint (x);
+%! %B = double (sinint (y));
+%! B = 1.02868521867373;
+%! assert (A, B, -5e-15);
+
+%!test
+%! %y = [2 3 sym(pi); exp(sym(1)) 5 6];
+%! x = [2, 3, pi; exp(1), 5, 6];
+%! A = sinint (x);
+%! %B = double (sinint (y));
+%! B = [1.60541297680269, 1.84865252799947, 1.85193705198247e+00; ...
+%!      1.82104026914757, 1.54993124494467, 1.42468755128051e+00];
+%! assert (A, B, -5e-15);
+
+%!assert (sinint (0), 0)
+%!assert (sinint (inf), pi/2)
+%!assert (sinint (-inf), -pi/2)
+
+%%tests against maple
+%!assert (sinint (1), 0.9460830703671830149414, -2*eps)
+%!assert (sinint (-1), -0.9460830703671830149414, -2*eps)
+%!assert (sinint (pi), 1.851937051982466170361, -2*eps)
+%!assert (sinint (-pi), -1.851937051982466170361, -2*eps)
+%!assert (sinint (300), 1.5708810882137495193, -2*eps)
+%!assert (sinint (1e4), 1.5708915453859619157, -2*eps)
+%!assert (sinint (20i), 1.2807826332028294459e7*1i, -2*eps)
+
+%!test
+%! x = (0:4) ';
+%! y_ex = [0
+%!         0.946083070367183015
+%!         1.60541297680269485
+%!         1.84865252799946826
+%!         1.75820313894905306];
+%! assert (sinint (x), y_ex, -4*eps);
+
+%!test
+%! x = -(0:4) ';
+%! y_ex = - [0
+%!           0.946083070367183015
+%!           1.60541297680269485
+%!           1.84865252799946826
+%!           1.75820313894905306];
+%! assert (sinint (x), y_ex, -4*eps);
+
+%!test
+%! x = 1i * (0:4).';
+%! y_ex = [0
+%!         1.05725087537572851*I
+%!         2.50156743335497564*I
+%!         4.97344047585980680*I
+%!         9.81732691123303446*I];
+%! assert (sinint (x), y_ex, -4*eps);
+
+%!test
+%! x = - 1i * (0:4).';
+%! y_ex = - [0
+%!           1.05725087537572851*I
+%!           2.50156743335497564*I
+%!           4.97344047585980680*I
+%!           9.81732691123303446*I];
+%! assert (sinint (x), y_ex, -4*eps);
+
+%!test
+%! % maple:
+%! % > A := [1+2*I, -2 + 5*I, 100, 10*I, -1e-4 + 1e-6*I, -20 + I];
+%! % > for a in A do evalf(Si(a)) end do;
+%! x = [1+2i; -2+5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i];
+%! A = [ 1.6782404878293681180 + 2.0396845546022061045*1i
+%!      -18.154174221650281533 + 1.6146414539230479060*1i
+%!       1.5622254668890562934
+%!       1246.1144901994233444*1i
+%!      -0.000099999999944461111128 + 0.99999999833338888972e-6*1i
+%!      -1.5386156269726011209 - 0.053969388020443786229*1i ];
+%! B = sinint (x);
+%! assert (A, B, -6e-16)