2
$\begingroup$

I'm trying to solve a nonlinear system of algebraic (not differential) equations which involves double integrals. I want to use Matlab. There seems to be two functions that can do double integration: dblquad and quad2d.

I noticed that the latter is much faster than the former (100+ times). But I don't know the accuracy. I used both functions to evaluate my integrals (which involves bivariate normal density) and got different answers. I chose a rectangular region around the mean in both cases, since Matlab (I think most numerical algorithms) seems to have difficulty finding the substantial values of the integrand if the region is too big and gives me 0.

My question is: does anyone know the difference of these two functions (in terms of integrating mechanism) as well as accuracy?

Thanks.

  • 0
    Did you make sure you exploited any symmetry in your integrands? The problem with the cubature of products of Gaussians is that if the peak is too small, sampling might miss it, which is supposed to contribute to the bulk of the integral. Could you maybe post the function you're trying to integrate?2011-05-03
  • 0
    What I'm trying to integrate is $E[c((a_1X_1)^{\alpha_1}+(a_2X_2)^{\alpha_2})X_1^{\alpha_1}]$ where $(X_1,X_2)$ is bivariate normal (correlated) and $c$ is some piecewise linear/quadratic function.2011-05-03
  • 0
    Note that `dblquad` is essentially a Cartesian product version of `quad` (i.e., think of computing the double integral by using `quad` with an integrand whose expression also involves `quad`), and thus is likely to be slow (due to the adaptive nature of `quad`) if the integrand has a lot of "features".2011-05-03
  • 0
    Anything special about the $a_i$ and the $\alpha_i$?2011-05-03
  • 0
    I double checked using Mathematica function NIntegrate[]. In that case, I simply integrate over the whole plane and don't have to choose a smaller integration region. It seems that Mathematica agrees with dblquad. But that too slow for solving my system because it calls the function many times.2011-05-03
  • 0
    No, nothing special, other than $a_i>0$, \alpha_i>1$.2011-05-03
  • 0
    Yes, as I said, you can expect `dlbquad` to be slow since it slices your integrand, applies `quad` to each of those slices, and then uses `quad` again to properly bring together those slices. *Mathematica*'s `NIntegrate[]` for multidimensional integrals uses a different cubature rule from MATLAB; if `dblquad`'s results agree with it, then you can have some confidence in the correctness. I don't know what algorithm `quad2d` uses; let me get back to you on that.2011-05-03
  • 0
    If the $\alpha_i$ are nonintegral, I'd expect the integrand to have a few singularities that may be affecting the results you're getting.2011-05-03
  • 0
    $\alpha_i$ are not integers, and I only integrate over $[0,2]\times[0,2]$ because the means are 1,1, and standard deviations are only around 0.12011-05-03
  • 0
    I went back and added an additional parameter ('AbsTol',1e-10) in quad2d and get an answer reasonably close to dblquad. I couldn't find the default AbsTol but I guess that's because of low accuracy requirement by default. With that higher accuracy, it's still very fast.2011-05-03
  • 0
    Ah, so my guess that `quad2d` didn't sufficiently sample your function was correct.2011-05-03

0 Answers 0