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

4 次查看（过去 30 天）
David Zhang
David Zhang 2012-2-3
编辑： Matt J 2013-9-28
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.

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

采纳的回答

Walter Roberson
Walter Roberson 2012-2-3
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 个评论
显示 6更早的评论
David Zhang
David Zhang 2012-2-4
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
Walter Roberson 2012-2-4
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.

请先登录，再进行评论。

更多回答（1 个）

Mike Hosea
Mike Hosea 2012-2-6
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 个评论
显示 3更早的评论
David Zhang
David Zhang 2012-2-19
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
Mike Hosea 2012-2-23
Use QUADGK and here again, put the singularity on the boundary: a = quadgk(f,-2,0) + quadgk(f,0,2)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Special Values 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by