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!
3 个评论
Jan
2013-7-4
"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
2013-7-4
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
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 个评论
Vera
2013-7-8
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
2013-7-8
0 个投票
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
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
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
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
