I'm trying to numerically evaluate this integral in Matlab using quad and dblguad functions, but I'm running into problems due to the way the inner integrals are structured. I cant use dblquad (double integration) to integrate over $u$ and $v$ because the function in the innermost integral is a function of $u$, $v$ and $z$. I guess my question is there any way I can get $z$ out of the two inner integrals. This way I can evaluate the inner integrals and all is left is just a function of $z$ for the outermost integral.

This is my code, but it's not working, because of the problem I mentioned above.

F1 = @(u,v) u.*v.*exp(-0.5.*u.^2)./(1+2*z.*u.*v);  
F2 = @(v) v;
F3 = @(z) exp(-z)./sqrt(z);
I1 = dblquad(F1,0,1e5,2,1e5);
I2 = quad(F2,2,1e5);

First observe that

\begin{align} \int_a^\infty \int_0^\infty v e^{-u} - \frac{v u e^{-u^2/2}}{1 + 2 u v z} \, du \, dv &= \int_a^\infty v \int_0^\infty e^{-u} \, du - \int_0^\infty \frac{v u e^{-u^2/2}}{1 + 2 u v z} \, du \, dv \\ &= \int_a^\infty v \cdot 1 - \int_0^\infty \frac{v u e^{-u^2/2}}{1 + 2 u v z} \, du \, dv. \end{align}

This allows us to rewrite the inner two integrals as a standard double integral:

$$\int_0^{\infty}\frac{e^{-z}}{\sqrt{z}} \exp\left[-\int_a^\infty \left( v - \int_0^\infty \frac{vue^{-\frac{u^2}{2}}}{1+2uvz}\,du \right) \,dv\right]\,dz = \int_0^{\infty}\frac{e^{-z}}{\sqrt{z}} \exp\left[-\int_a^\infty \int_0^\infty v e^{-u} - \frac{vue^{-\frac{u^2}{2}}}{1+2uvz}\,du\,dv\right]\,dz.$$

Next, in MATLAB we define a function myfun.m which, given a vector of $z$-values, evaluates the function

$$z \mapsto \frac{e^{-z}}{\sqrt{z}} \exp\left[-\int_a^\infty \int_0^\infty v e^{-u} - \frac{vue^{-\frac{u^2}{2}}}{1+2uvz}\,du\,dv\right].$$

function f = myfun(z)

a = 1;
max_x = 1000;
max_y = max_x;

f = zeros(size(z));

for idx = 1:numel(z)
    f(idx) = integral2(@(u, v) v.*exp(-u) - v.*u.*exp(-u.^2/2)./(1 ...
                               + 2*u.*v*z(idx)),0,max_x,a,max_y);

f = exp(-z - f)./sqrt(z);

The reason to use max_x and max_y instead of inf is because otherwise the integral diverges and MATLAB throws an error. (This might be an indication that the problem is not well-posed; I'm not sure.) It might be possible to vectorize the loop, but it isn't clear to me how to do that.

Lastly, q = integral(@myfun, 0, inf) will compute the desired integral (1).

  • $\begingroup$ Thanks for the answer. Can I use quad instead of integral2 and integral, I'm using an older version of Matlab. Also myfun takes a vector of z values as input, where are these values created?. $\endgroup$ – user4259 Oct 31 '13 at 5:07
  • $\begingroup$ Glad to help! You could replace integral with quadgk and integral2 with dblquad (no change is necessary to the values inputted into the functions). Let me know if this works with your version of Matlab. (The reason for quadgk instead of quad is due to the improper integral - quadgk can compute improper integrals whereas quad cannot.) The vector of z values are created by the integral function (or quadgk function) when it calls @myfun. $\endgroup$ – Kyle Oct 31 '13 at 6:04
  • $\begingroup$ Yes it works like a charm, thank you. However, I'm having trouble choosing the proper values of the upper limits in dblquad. I've tried different values and I'm getting different answers for each, so I don't know where does it converge. Is there any alternative to dblquad that would accept inf as an upper limit directly like quadgk. $\endgroup$ – user4259 Oct 31 '13 at 22:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.