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
Best Answer