why dblquad can not be used to evaluate the bessel function of the second kind bessely

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.

Answer 1

[Walter Roberson]

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.

[David Zhang]

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]

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]

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]

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]

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]

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]

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]

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.

Answer 2

[Mike Hosea]

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

[David Zhang]

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]

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]

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]

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]

Use QUADGK and here again, put the singularity on the boundary: a = quadgk(f,-2,0) + quadgk(f,0,2)