How do I properly write this integration in Matlab?

Question

I'm trying to evaluate the following integral in matlab: http://i.imgur.com/Iuc4VT5.png

Here is my code:

alpha = 2;
F1 = @(u,v) 2*u.*v.*exp(-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);
quad(F3*exp(-(I2-I1)),0,1e5);

I'm getting the errors shown below. These errors don't show much, but I'm guessing it's because of the way I wrote F1. I defined F1 as a function of u and v for the double integral, but there is also a variable z which is the variable of the outer integral. I did that because there is no way I can separate z from the inner integrals. Is there any better way to write this integration?

Error in ==> @(u,v)2*u.*v.*exp(-u.^2)./(1+2*z.*u.*v)


Error in ==> dblquad>innerintegral at 73
fcl = intfcn(xmin, y(1), varargin{:}); %evaluate only to get the class below

Error in ==> quad at 76
y = f(x, varargin{:});

Error in ==> dblquad at 53
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

I'm choosing 1e5 to represent infinity.

You are telling us where the error is raised, but not what the error message is. — Daniel
@DanielR These are all the error messages I'm seeing in Matlab. — Thomas Jr
Well, running it on my system there surely is another message: Undefined function or variable 'z'. — A. Donda
@A.Donda I don't get that error, I'm running an older version of Matlab. But that could possibly be the error, I believe the way I defined F1 is not correct, but I don't know of any other way. Any suggestions? — Thomas Jr
@A.Donda Thanks. I'll see if I can get any help from math.stackexchange.com — Thomas Jr

A. Donda A. Donda · Accepted Answer · 2013-10-29T19:58:36

After having successfully answered a follow-up question by the same poster, I realize a relevant part of this answer is wrong. I'd delete the answer, but I can't, since it is accepted. Therefore this disclaimer...

The simple answer is: Your definition of F1 contains a reference to z, but it is not specified as an argument of that function.

However, it won't help to specify z as an additional argument, because then I1 is no longer a constant, but itself a function of z.

I'm not an expert on numerical integration, but as far as I can see this means that you cannot numerically integrate your expression, at least not using a combination of quad and dblquad. The argument to the outer exponential function is just not a constant, and a numerical integration cannot return a function.

It might be possible to rearrange the integral to bring it into a form that can be numerically integrated, but I can't tell you how.

Another problem is that "representing" infinity by 10^5 is not necessarily a useful approximation – it all depends on the behavior of the function that is being integrated. A possible trick might be to do a variable substitution such that each variable that goes to infinity is written as a function of another variable with a limited range.

My recommendation: Try to get as far as possible evaluating this integral analytically, and only use numerics when you are sure there is no analytical approach. And try to get help on that on math.stackexchange.com, because it is not a programming problem.

How do I properly write this integration in Matlab?

1 Answers