why dblquad can not be used to evaluate the bessel function of the second kind bessely
David Zhang asked
on 3 Feb 2012
Latest activity:
Answer by Mike Hosea
on 6 Feb 2012
Hi All, I like to calculate the double integral of Bessel function as clear clc a=0; b=2; syms x y A= dblquad(@(x, y) -abs(x-y).*bessely(1, abs(x-y)), a, b, a, b); where the error is noted as Warning: Infinite or Not-a-Number function value encountered. when i take the quad of Bessel function as clear clc a=0; b=2; syms x A= quad(@(x, y) -abs(x).*bessely(1, abs(x)), a, b); The results is OK! Anyone know what is the problem with my dblquad? Thank you a lot. Products |
|---|
2 answers
Walter Roberson answered
on 3 Feb 2012
Accepted answer
You do not need the "syms" line at all.
Your integration ranges are the same for x and y, so there are places where x == y and at those places, abs(x-y) will be abs(0).
BesselY(1,x) is
sum( (1/2) * (-2*Psi(2+_k1) + 1/(_k1+1) + 2*ln(x) - 2*ln(2))*x^(1+2*_k1)*2^(-2*_k1)*(-1)^(3*_k1) / (factorial(_k1+1)*factorial(_k1)*Pi), _k1 = 0 .. infinity) - 2/(x*Pi)
and if we examine even just the last term of that, 2/(x*Pi) we can see that BesselY(1,x) is not going to be a finite number.
I do see you multiplying the bessely term by abs(x-y) but if the bessely term has already generated inf or nan, multiplying that by 0 is not going to give you 0.
Why do you not get it in the simpler case, considering that one of your bounds is 0? I do not know for sure -- but in the longer case, there are multiple (many!) points where you would evaluate at 0, not just a single one.
8 comments
Show
5 older comments
David Zhang
on 3 Feb 2012
Hi Walter, thanks for your help. I do change the bounds as not 0, but the same error is also noted. can you give me a simple note of computing the double integration of Bessel function?
Thanks
Walter Roberson
on 3 Feb 2012
Single integral of your first function turns out to have an analytic expression, but I do not know yet if the double integral does. I am running it through some probe routines now.
David Zhang
on 3 Feb 2012
Yes, you are correct, for the single integration, it is
int(abs(x)*BesselY(1, abs(x)), x = 0 .. 1)=(1/2)*Pi*(BesselY(1, 1)*StruveH(0, 1)-BesselY(0, 1)*StruveH(1, 1))
I calculate that in Maple, but for the double integration, i fail to get the result.
I am wild of how to get the double integration numerically. Thanks
Walter Roberson
on 3 Feb 2012
MuPAD does not appear to know the StruveH function, so in Maple, convert(%,hypergeom)
Then, take the expression and wrap it around with a quad() in order to do the numeric integral over the second variable.
David Zhang
on 3 Feb 2012
Walter, Thanks, but i have to do it using matlab. It is just not a small case to know the result of integral. Do you know how to get that using Matlab? Thanks again.
Walter Roberson
on 3 Feb 2012
a=0;
b=2;
syms x y
intx = matlabFunction( int(-abs(x-y).*bessely(1, abs(x-y)), x, a, b), y);
dblint = quad(intx, a, b);
David Zhang
on 4 Feb 2012
Walter, the code is tested with the error and warning
Warning: Explicit integral could not be found.
> In sym.int at 64
??? Error: The expression to the left of the equals sign is not a valid target for an assignment.
BTW, intx = matlabFunction( int(-abs(x-y).*bessely(1, abs(x-y)), x, a, b), y) seems to take so long time.
Walter Roberson
on 4 Feb 2012
I do not have the symbolic toolbox to test with (just Maple); it could be that it is too complex for MuPAD to work with.
Unfortunately MuPAD lacks the convert() operation; the closest it has is coerce() which has a different intent.
If the first level integration must be performed within MATLAB, and the MuPAD int() is not able to do it, I do not see what can be done other than numeric integration.
Perhaps if you defined
function r = valhere(x,y) %x can be vector, y is scalar
r = zeros(size(x));
idx = abs(x-y) <= 1e-6; %use appropriate tolerance
r(idx) = -abs(x(idx) - y) .* bessely(1, abs(x(idx)-y));
end
dblint = dblquad(@valhere, a, b, a, b);
If you do not know ahead of time what the interior of the function f(x,y) looks like,
function r = evalhere(f,x,y)
r = f(x,y);
r(isnan(r)) = 0;
r(isinf(r)) = 0; %a bit dubious really
end
dblint = dblquad(@(x,y) evalhere(f,x,y), a, b, a, b);
However, if your function is known in advance, such as the specific function here, then I would say integrate it once in Maple, convert() it to hypergeom, copy-and-paste into MATLAB changing the BesselY to bessely and the like, then quad() that.
Mike Hosea answered
on 6 Feb 2012
Integrate using QUAD2D and put the singularity on the boundary.
>> f = @(x,y)-abs(x-y).*bessely(1,abs(x-y))
f =
@(x,y)-abs(x-y).*bessely(1,abs(x-y))
>> A = quad2d(f,0,2,@(x)x,2) + quad2d(f,0,2,0,@(x)x)
A =
2.81899103740708
5 comments
Show
2 older comments
David Zhang
on 6 Feb 2012
Hi Mike, it is done. Thanks a lot. BTW, the matlab manual said that quad2d can only deal with the 2-d plane problems, my question is: if the integral are multidimensional, how this problems can be solved then?
thanks!
Mike Hosea
on 7 Feb 2012
QUAD2D integrates over "planar regions" that are defined by a <= x <= b and c(x) <= y <= d(x). The planar region is the "support". The integral that is performed is a double integral. DBLQUAD can only integrate over rectangular regions a <= x <= b and c <= y <= d. QUAD2D is better than DBLQUAD on most problems.
Mike Hosea
on 9 Feb 2012
Maybe I misunderstood your question. If you mean, "How do I solve problems like this in in higher dimensions?" then the only answer I have for you is to nest QUAD2D and/or QUADGK using the arrayfun function. This is difficult, but I believe I have supplied codes in my other answers to evaluate triple integrals using QUADGK with QUAD2D, and maybe I have even supplied code to evaluate a 4-D integral nesting QUAD2D with QUAD2D. Clearly this will get very expensive very fast and may not be practical in more than a 3 or 4 dimensions.
David Zhang
on 19 Feb 2012
Hi Mike,
I have another question to ask you, abt the bessely function. Simply for the 1-D integration with
f=@(x) -abs(x).*bessely(1, x); %%(***)
a=quad(f, -2, 2)
This code my induce the inf or NAN because of bessely(1, 0) are infinite number, this can be resolved using function quad1d as you previously said. My new problem is to numerically integrate (***) by
using Gaussian quadrature, so when the quadrature point are 0, so how i can solve it now? because the limit{x-->0} of bessely(1, x) are not existed. Thank you!
Mike Hosea
on 23 Feb 2012
Use QUADGK and here again, put the singularity on the boundary: a = quadgk(f,-2,0) + quadgk(f,0,2)
Contact us at files@mathworks.com
0 comments