Double integral over an anonymous function

I am trying to do the following double integral:

where

I attempt this by

q = 0.01:0.01:50; % range is somewhat arbitrary
fun = @(q) integral2(@(u,sigma) (((sin(q.*(B/2)*sqrt(1-sigma.^2).*cos(pi*u./2))./(q.*(B/2)*sqrt(1-sigma.^2).*cos(pi*u./2))).*(sin(q.*(B/2)*sqrt(1-sigma.^2).*a*sin(pi*u./2))./(q.*(B/2)*sqrt(1-sigma.^2).*a*sin(pi*u./2)))).^2).*(sin(q.*B*c*sigma./2)./(q.*B*c*sigma./2)).^2,0,1,0,1);
pq = arrayfun(@(q)fun(q),q);

I imagine I made a mistake, especially regarding the element wise operations, but I don't know how to write this in a cleaner way because that breaks up the anonymous function. Can somebody help me clean this up?

Answer 1

[Walter Roberson]

We could probably produce a more compact function if we knew the sign() of the various constants

syms a B c mu q sigma x real
assume(a ~= 0)
assume(B ~= 0)
assume(c ~= 0)
assume(q > 0)     %an important restriction that avoids division by 0
assume(sigma ~= 0)

Pi = sym(pi);
S(x) = sin(x)/x

S(x) = 

inner_phi = (S(mu/2*cos(Pi/2*mu)) * S(mu*a/2 * sin(Pi/2*mu)))

inner_phi = 

phi(mu,a) = int(inner_phi^2, mu, 0, 1)

phi(mu, a) = 

Pmu_inner = phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2

Pmu_inner = 

Pmu(mu) = int(Pmu_inner, sigma, 0, 1)

Pmu(mu) = 

P(q) = Pmu(q*B)

P(q) = 

Pfun = matlabFunction(P, 'vars', [q, a, B, c, sigma])

Pfun = function_handle with value:
    @(q,a,B,c,sigma)integral(@(sigma)1.0./B.^2.*1.0./c.^2.*1.0./q.^2.*1.0./sigma.^2.*sin((B.*c.*q.*sigma)./2.0).^2.*integral(@(mu)1.0./a.^2.*1.0./mu.^4.*sin((a.*mu.*sin((mu.*pi)./2.0))./2.0).^2.*1.0./cos((mu.*pi)./2.0).^2.*1.0./sin((mu.*pi)./2.0).^2.*sin((mu.*cos((mu.*pi)./2.0))./2.0).^2.*1.6e+1,0.0,1.0).*4.0,0.0,1.0)

Is this going to be usable in practice?

Possibly not. You have S() of expressions involving sin(mu) where mu is q*B and q ranges 0 to 50. Unless B is in the range 0 to 1/50 or so, then you are going to be taking sin() of an expression that ranges further than pi/2 so you are going to (at least in theory) be getting a 0 from that part. But it is an argument to S which is sin(x)/x and the x part is the sin(q*B*pi/2) and if that inner sin() can go to 0, then you are going to have division by 0, which is going to give you problems during numeric evaluation.

The sinc() function itself does not have problems with division by 0 because mathematically it is not just plain sin(x)/x but rather limit(sin(x+delta))/(x+delta) as delta->0 and the limit() goes to 1 there because near 0 sin(x) can be closely approximated by x, and limit (x+delta)/(x+delta) goes to 1.