Implementing generalized cosine and sine integrals

      hello everybody,
      I have been trying to implement a formula which includes generalized sine and cosine integrals. I didn't succeed yet. 
      Any of you please might help me with this??
      The formula is: 

       W = 2[ sinh^-1 (h/a) - C(2ba, 2bh) - jS(2ba, 2bh) ] 
       + (j/bh) *(1 - e^(-jbh)),
      where b = 57.8, h = 25 mm (0.025 m), a = 0.5 mm (0.0005 m).

      C(x, y) and S(x, y) are the generalized cosine and sine integrals defined as:

      C(x, y) = Integral [ cos y / t^y  dt], with 0-x as integration borders
      S(x, y) = Integral [ sin y / t^y  dt], with 0-x as integration borders

    I have been trying to calculate the S(x, y) like (same for the C(x, y)):
    b = 57.8;
    h = 0.025; % 25 mm
    a = 0.0005; % 5 mm

    k = 2* b * h;
    l = 2 * b * a;
    fun = @(x) sin(x) ./ (x.^k);
    q = integral (fun, 0, l);

but it doesn't work.

Can, please, any of you help me with this?
Thanks in advance for your time!

[Jan]

"It doesn't work" does not allow to reconsider which problem you are faced with. Do the results differ from your expectations? Do you get an error message and if so, which one? Do the error message reveal important information, e.g. that integral operates on symbolic variables only?

[Vera]

    hi Jan,
    thanks for answering.
    This is what it says:
    Warning: Infinite or Not-a-Number value encountered. 
  > In funfun\private\integralCalc>iterateScalarValued at 349
    In funfun\private\integralCalc>vadapt at 133
    In funfun\private\integralCalc at 76
    In integral at 89
    In integral_gener_sin at 8 

Would you please tell me what is the problem? I am not sure I understand it...
Thank you very much!

[Oliver K]

It was just the warning in the output. Didn't you get any result in q? What was your expected result of q?

Answer 1

[Oliver K]

The approach suggested by Mike is a good one. I modified Mike's functions a little to accommodate the case where k >= 2

%%Generalized integral of cos(x)/x^k
function q = gencosint(b,k,varargin)
  if k < 1
      q = quadsas(@cos,0,b,-k,0,varargin{:});
  elseif k == 1
      q = integral(@(x) cos(x)./x,0,b,varargin{:});
  else
      q = -(gensinint(b,k-1,varargin{:}) + cos(b)/(b^(k-1)))/(k-1);
end

%%Generalized integral of sin(x)/x^k 
function q = gensinint(b,k,varargin)
  if k < 1
      q = quadsas(@sin,0,b,-k,0,varargin{:});
  elseif k == 1
      q = integral(@(x)sin(x)./x,0,b,varargin{:});
  else
      q = (gencosint(b,k-1,varargin{:}) - sin(b)/(b^(k-1)))/(k-1);
end

%%example from the question
b = 57.8;
h = 0.025; % 25 mm
a = 0.0005; % 5 mm
k = 2* b * h;
l = 2 * b * a;
q = gensinint(l,k)

[Vera]

Thank you Oliver!

[Oliver K]

You are welcome. I think that I made a mistake in previous answer with the integration by parts. When doing the integration by part for cos(x)/x^k, evaluating the value of cos(x)/x^k with x=0 we would have 1/0=Inf. Similarly, for sin(x)/x^k we have 0/0.

You may want to double-check your input values though. As far as I can see, your input for k is 2.8900. Generalized integral of sine and cosine don't converge in this case.

Answer 2

[Walter Roberson]

At 0, your formula comes out to 0 / 0 which is Not A Number.

integral() is supposed to be able to handle a singularity at the lower limit, so it should just be a warning.

If your results are not good you might want to look on the integral() documentation page to see how to specify tolerances.

Answer 3

[Mike Hosea]

Well, first of all, I think we need 0 < k < 2 for the generalized sine integral and 0 < k < 1 for the generalized cosine integral. But the main issue is that the singularity is too strong. Here is what I would suggest. Get QUADSAS

http://www.mathworks.com/matlabcentral/fileexchange/32542-quadsas

This should evaluate the generalized cosine integrals directly as well as the generalized sine integrals where k < 1. For 1 < k < 2 you can employ integration by parts and reduce the generalized sine integral to a generalized cosine integral. The "varargin" parts here are so you can pass 'AbsTol' and 'RelTol'.

function q = gencosint(b,k,varargin)
q = quadsas(@cos,0,b,-k,0,varargin{:});

function q = gensinint(b,k,varargin)
if k < 1
    q = quadsas(@sin,0,b,-k,0,varargin{:});
elseif k == 1
    q = integral(@(x)sin(x)./x,0,b,varargin{:});
else
    q = (gencosint(b,k-1,varargin{:}) - sin(b)/(b^(k-1)))/(k-1);
end

[Vera]

Thank you very much for your time. Very helpful.