how to improve accuracy of double integral (i.e. quad2d, iterated integral-quadgk twice) when variables are close to zero?

given rho,alpha,theta0,r0 are fixed constant, I want to do double integral of a complicated function f(s,t) from (0,inf,0,inf). f(s,t) is listed as follows:

%when 0<s<t

term1 = @(s,t) pi*sin(alpha)/2/alpha^2./sqrt(s.*(t-s*cos(alpha)^2))./(t-s);

term2 = @(s,t) exp(-abs(real(r0^2./2.*(1-cos(2.*alpha))./((t-s)+(t-s.*cos(2.*alpha))))));

term3 = @(n,s,t) n.*sin(n.*pi*(alpha-theta0)/alpha)... .*besseli(n.*pi/2/alpha,r0^2./2./s.*(t-s)./((t-s)+(t-s.*cos(2*alpha))),1); term3_sum = @(s,t) arrayfun(@(si,ti) sum(term3(maxmode_1,si,ti)),s,t);

f = @(s,t) term1(s,t).*term2(s,t).*term3_sum(s,t)./s./t;

% when s>t>0

term1_alt = @(s,t) pi*sin(alpha)/2/alpha^2./sqrt(t.*(s-t*cos(alpha)^2))./(s-t);

term2_alt = @(s,t) exp(-abs(real(r0^2./2.*(1-cos(2.*alpha))./((s-t)+(s-t.*cos(2.*alpha))))));

term3_alt = @(n,s,t) n.*sin(n.*pi.*theta0/alpha)... .*besseli(n.*pi/2/alpha,r0^2./2./t.*(s-t)./((s-t)+(s-t.*cos(2*alpha))),1);

term3_sum_alt =@(s,t)arrayfun(@(sii,tii)sum(term3_alt(maxmode_1,sii,tii)),s,t);

f= @(s,t) term1_alt(s,t).*term2_alt(s,t).*term3_sum_alt(s,t)./s./t;

I know that I have singularity along the line s=t. I break the double integral from lower triangle and upper triangle to estimate the total double integral. I used both two methods to do double integral i.e. quad2d and iterated integral. I used iterated integral(quadgk twice) method to double integral from (0,inf,s,inf)+(0,inf,0,s). However, it fails to integral from 0 since s,t close to zero, f will converge to inf. I only can put lower bound equal to 10^-3. HOwever,both methods gives me solution of 0.0043 while I am expecting my true value is 0.0044. I think my error comes from the double integral area where s,t are close to zero. When s,t are close zero, my function converges to inf (since ./s./t). quad2d could lose accuracy when s,t close to zero. Iterated integral method fails to integral from 0. Do any one know how to better estimate the double integral of f when s,t are around zero (f will converge to inf)? How to improve my accuracy of quad2d or iterated integral method when s,t are around zero?

Thank you in advance,

Regards,

Lu

Answer 1

[Mike Hosea]

Perhaps the difference here is just that you are missing the area in the square [0,.001]x[0,.001]. It would be easier to look into this if you would provide a working example with values for all parameters. However, if the singularity is too strong for these numerical methods, it won't matter. A common approach to deal with that situation is to develop a series approximation in a neighborhood of the singularity. You use the series approximation over the small neighborhood where it can give an accurate result. Then you can use the numerical method from the edges of that neighborhood, as you have already done.

Answer 2

[Lu]

Thank you, Mike. It really helps.

Sincerely,

Lu