"kumar vishwajeet" <firstname.lastname@example.org> wrote in message news:email@example.com... > I am integrating the following function between L and U using dblquad. I > get "NaN" and a warning of singularity. F = > real(constantForM.*((x-L).^(A(1)-1).*(y-L).^(A(2)-1).*(U-L-((x-L)+(y-L))).^(A(3)-1)).^2 > > where, constantForM = 3.378137617443966e-038. > L = 1e-7 > U = 3.2e5 > A = [3.75 0.25 0.02] > > In order to check for singularity, I generated 1e6 points between L and U > and evaluated F at each of the points. I get a sharp peak(singularity) at > halfway. But the value at that point is still less than 1. In fact the > maximum value of F along these points is 0.00116 and minimum is close to > zero. Then why do I get NaN??
When y = L (assuming you're integrating over the square [L, U] x [L, U] the middle term in your integrand is 0^(A(2)-1) = 0^(-0.75). That's Inf. If x is also equal to L, the first term in the product is 0^(2.75) = 0. Then 0*Inf = NaN and once one NaN value creeps into a calculation, it propagates.
> Another thing: > I used Monte Carlo Integration using the following method: > 1. Evaluated F at all 1e6 points. > 2. Found the average of all those values. > 3. Multiplied by the volume i.e. (U-L)*(U-L). > 4. I get 119 as the answer. > > Which of these two methods is correct??
Your integrand isn't defined in the lower-left corner of your region of integration, if my assumption about your region of integration is correct. Therefore I'd say the NaN result you received is correct.