Implementing generalized cosine and sine integrals

14 次查看（过去 30 天）
Vera
Vera 2013-7-4
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!
  3 个评论
显示 1更早的评论
Vera
Vera 2013-7-4
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
Oliver K 2013-7-7
It was just the warning in the output. Didn't you get any result in q? What was your expected result of q?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Oliver K
Oliver K 2013-7-8
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)
  2 个评论
Oliver K
Oliver K 2013-7-11
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.

请先登录，再进行评论。

更多回答（2 个）

Walter Roberson
Walter Roberson 2013-7-8
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.
Mike Hosea
Mike Hosea 2013-7-8
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

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by