Mercurial > octave
view scripts/specfun/sinint.m @ 27919:1891570abac8
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2020.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Mon, 06 Jan 2020 22:29:51 -0500 |
parents | b442ec6dda5c |
children | bd51beb6205e |
line wrap: on
line source
## Copyright (C) 2017-2020 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this distribution ## or <https://octave.org/COPYRIGHT.html/>. ## ## ## 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 ## <https://www.gnu.org/licenses/>. ## -*- 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 if (! isnumeric (x)) error ("sinint: X must be numeric"); endif ## Convert to floating point if necessary if (isinteger (x)) x = double (x); endif ## Convert to column vector orig_sz = size (x); x = x(:); if (iscomplex (x)) ## Work around reshape which narrows to real (bug #52953) x = complex (real (x)(:), imag (x)(:)); else x = x(:); end ## Initialize the result y = zeros (size (x), class (x)); tol = eps (class (x)); todo = true (size (x)); ## Special values y(x == 0) = x(x == 0); # correctly signed zero y(x == Inf) = pi / 2; y(x == - Inf) = - pi / 2; todo = ((todo) & (x != 0) & (x != Inf) & (x != - Inf)); ## For values large in modulus we use the relation with expint flag_large = abs (x) > 2; xx = x(flag_large & todo); 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(todo & flag_large) = yy; ## For values small in modulus we use the series expansion todo = (todo) & (! flag_large); xx = x(todo); ssum = xx; # First term of the series expansion yy = ssum; flag_sum = true (nnz (todo), 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) >= tol); it++; endwhile y(todo) = yy; y = reshape (y, orig_sz); endfunction %!assert (sinint (1.1), 1.02868521867373, -5e-15) %!test %! x = [2, 3, pi; exp(1), 5, 6]; %! A = sinint (x); %! B = [1.60541297680269, 1.84865252799947, 1.85193705198247e+00; ... %! 1.82104026914757, 1.54993124494467, 1.42468755128051e+00]; %! assert (A, B, -5e-15); ## Test exceptional values %!assert (sinint (0), 0) %!assert (signbit (sinint (-0))) %!assert (sinint (Inf), pi/2) %!assert (sinint (-Inf), -pi/2) %!assert (isnan (sinint (NaN))) ## Check single data type is preserved %!assert (class (sinint (single (1))), "single") ## Tests against Maple %!assert (sinint (1) , 0.9460830703671830149414, -2*eps) %!assert (sinint (-1) , -0.9460830703671830149414, -2*eps) %!assert (sinint (pi) , 1.851937051982466170361, -3*eps) %!assert (sinint (-pi), -1.851937051982466170361, -3*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, -3*eps) %! B = sinint (single (x)); %! assert (A, B, -3*eps ("single")) ## FIXME: Need a test for bug #52953 %#!test <*52953> ## Test input validation %!error sinint () %!error sinint (1,2) %!error <X must be numeric> sinint ("1")