[Scilab-users] Si function

[Federico Miyara]

Dear All,

Just in case somebody is interested, find attached a Scilab function to 
compute the sine integral function Si (the integral from 0 to x of the 
sinc function), which cannot be expressed in closed form with elementary 
functions.

/It is preliminary, it doesn't test for appropriate input argument.

It works for real or complex matrices or N-D arrays.

It can be easily modified to get the cosine integral function.

Regards,

Fderico Miyara
/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20200206/bd0a5a4a/attachment.htm>
-------------- next part --------------
function y = Si(x)
    // This function computes the sine integral function, i.e., 
    // the integral from 0 to x of the sinc function sin(x)/x
    // Usage:
    //       y = Si(x)
    // where
    //       x: Any real or complex scalar, vector, matrix  
    //          or N-D aray
    //       y: The result of applying the sine integral to
    //          each component. Its size is equal to size(x)
    //
    // Reference: 
    // Abramowitz, M. and Stegun, I. (eds.) "Handbook of 
    // Mathematical Functions". Dover. New York, 1972. 
    // ---------------------------
    // Author: Federico Miyara
    // Date:   2020-02-06
        
    siz = size(x);
    
    // Indexes to values of the argument with absolute 
    // value not greater than 20
    I = find(abs(x)<=20);
    // Indexes to values of the argument greater than 20
    // and non-negative
    J = find(abs(x)>20 & x<>-abs(x));
    // Indexes to values of the argument smaller than 20
    // and negative 
    K = find(abs(x)>20 & x==-abs(x));
    
    // Apply algorithm based on Taylor series 
    y_small = Si_small(x(I));
    
    // Apply algorithm based on asymptotic approximation
    // This approximation is not valid for negative arguments,
    // but as Si is an odd function this is easy to fix 
    y_large_pos = Si_large(x(J));
    y_large_neg = -Si_large(-x(K));
    
    // Put all together
    y = zeros(x);
    y(I) = y_small;
    y(J) = y_large_pos;
    y(K) = y_large_neg;
       
endfunction

function y = Si_small(x)
    N = 30
    coe = (-1).^(0:N)./(1:2:2*N+1)./factorial(1:2:2*N+1); 
    Sipoly = poly(coe, "u", "c");
    y = x.*horner(Sipoly, x.^2);  
endfunction

function y = Si_large(x)
    M = 15
    // Coefficients of the Laurent series of the auxilary
    // functions for the asymptotic approximation 
    sig = (-1).^(0:M);
    cf = factorial(0:2:2*M).*sig;
    cg = factorial(1:2:2*M+1).*sig; 
    
    // Polynomials with the preceding coefficients
    pf = poly(cf,"u","c");   
    pg = poly(cg,"u","c");  

    // Computation of the auxiliary functions
    f = horner(pf, 1./x.^2)./x
    g = horner(pg, 1./x.^2)./x.^2
    
    // Computation of tha asymptotic approximation
    // fo large arguments
    y = %pi/2 - f.*cos(x) - g.*sin(x)
endfunction