How to do numerical multiple integral (more than triple) by using matlab?

Dear friends, my question is How to do numerical multiple integral (more than triple) by using matlab?

For example, I have a function,Y = Func. It has 5 var. Hence, Y= func(a,b,c,d,e). If I want to do the integration with respect to a,b,c,d,e. (var a is the first one, e is the last one.) And region of a = [0 ,b] (this is a function handle), b = [0,c], c= [0.d], d= [0,e], e=[0, Inf].

Thank you so much for your time.

Thank you.

[Mike Hosea]

You have failed to supply any arguments to out1fcn. Remember that function handles always require arguments when you are evaluating them.

Answer 1

[Mike Hosea]

You can use INT if your problem can be handled symbolically. If not, numerical integration of a 5-fold integral in MATLAB requires nesting INTEGRAL, INTEGRAL2, and INTEGRAL3. I got tired of explaining how to do this, so I wrote integralN and put it on the file exchange.

http://www.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

There are some examples in the text at the top of the file. Basically, it's a straightforward extension of INTEGRAL2 and INTEGRAL3 in terms of how you call it.

You're going to have to reorder your variables because of the way your region is defined. By convention, the outermost integral is the first variable, and the innermost integral is the last. You need that reversed.

q = integeralN(@(e,d,c,b,a)fun(a,b,c,d,e),0,inf,0,@(e)e,0,@(e,d)d,0,@(e,d,c)c,0,@(e,d,c,b)b)

[sun]

Dear Mike. Thank you so much. Please allow me to take a look. I will re-edit my question again soon. You have a nice day.

[Mike Hosea]

This runs. I have no idea when it will stop running, though. :) Notice that I chose to stay in the symbolic world until the last moment when I used matlabFunction to create the numerical integrand. That was so that I could try symbolic integrations to reduce the dimensionality. None of them worked for me. Note that to do this, you need 4 arguments to INT, and you don't need any function handles. For example if you wanted to integrate x^2 + y*x + 3 with respect to x and do it from x=0 to x=y, just do this

>> syms x y
>> q = int(x^2 + y*x + 3,x,0,y)
q =
(5*y^3)/6 + 3*y

----------------------------------------------------------

%%%==== just some parameters ====
a=4;
la1=1/(pi*500^2); la2= la1*5;
p1=25; p2=p1/25;
sgma2=10^(-11);
index=1;
g=2./a;
syms r u1 u2 u3 u4 u5
index = -2;
powe= index ;
seta= 10^powe;
xNor = ( (u5./u1).^(a./2)+ (u5./u2).^(a./2) + (u5./u3).^(a./2)+ (u5./u4).^(a./2) + 1    ).^(2./a);
x = (xNor).^(0.5) * seta^(-1/a);
q=pi.*(la1.*p1.^(2./a)+la2.*p2.^(2./a));
%%%====  parameters end ====
fun1 =  r./(1+ r.^a );
out1 = int(fun1, x, Inf) ;
ysym  = exp(-u3*(1+2*(out1)/((( (u5/u1)^(a/2)+ (u5/u2)^(a/2) + (u5/u3)^(a/2) + ...
    (u5/u4)^(a/2) + 1)^(2/a))*seta^(-2/a))))*exp(-sgma2*q^(-a/2)*seta*u3^(a/2)/ ...
    ((((u5/u1)^(a/2)+(u5/u2)^(a/2) + (u5/u3)^(a/2)+ (u5/u4)^(a/2) + 1)^(2/a))^(a/2)));
% You can try symbolically integrating ysym with respect to any u variable.
% My attempts were unsuccessful.
y = matlabFunction(ysym);
q = integralN(@(u5,u4,u3,u2,u1)y(u1,u2,u3,u4,u5),0,inf,0,@(u5)u5,0,@(u5,u4)u4,0,@(u5,u4,u3)u3,0,@(u5,u4,u3,u2)u2)

[Arny]

Hi Mike, Would it be possible to suggest how the following multidimensional integral could be coded in Matlab, so as to useintegralN.m?

integrate {p1, 0, a}, {p3, 0, b}, {p4, 0 c} func1 * func2 * func3 * func4 * func5

where

func1=((p3 p4)*(Sin (p1 r)))/p1

func2= p1p3 + p1p4 + p3p4

func3= (1/(p3^2 +1)^2)*(p4^2 +1)^2

func4=Log((p3+p1)^2/(p3-p1)^2)

func5=Log((p4+p1)^2/(p4-p1)^2)

Assumptions p1>0, p3>0, p4>0.

Kindly note that there are singularities that have to be avoided at p4=p1, p3=p1 etc.

Answer 2

[Shoaibur Rahman]

You can eliminate the variables one by one until you feel comfortable. You can make down to 3 or 2 or 1 variable before final integration, since Matlab has integral3, integral2 and integral to perform triple, double and single integration respectively.

y = @(a,b,c,d,e) a+b+c+d+e; % your function
y1 = @(a,b,c,d) integral(@(e) y(a,b,c,d,e),L,H); % L,H are the integration limits for e
y2 = @(a,b,c) integral(@(d) y1(a,b,c,d),L,H); % L,H are the integration limits for d
% and so on

Also look at integral doc for further input arguments if required.

[sun]

Dear sir, I think intergral is only for real or complex number. I can not give 'syms' for L, H.

[Shoaibur Rahman]

Try with this, although the explicit integral could not be found:

y1 = int(y, u1, 0, u2)
y2 = int(y1, u2, 0, u3)
y3 = int(y2, u3, 0, u4)
y4 = int(y3, u4, 0, inf)

Answer 3

[John D'Errico]

You need to consider the curse of dimensionality. Adaptive numerical integrations often take hundreds of function evaluations. Depending on the complexity of your functions, a single adaptive numerical integration might take hundreds or more function evaluations.

For example, a quick test of the integral function shows to integrate a simple function like sin(x) over the interval [0, 10*pi], took 390 function evaluations.

Even a very simple and smooth function like exp, over the interval [0, 2] took 150 function evaluations.

So suppose these were a reasonable counts for integral? A 5 level iterated integral could then take on the order of

(390)^5
ans =
   9.0224e+12

or

(150)^5
ans =
   7.5938e+10

function evaluations! In either case, we might be talking about something between 75 billion and 9 trillion function evaluations.

So if you have no choice but to do the multiple integration using such a numerical tool, do what you must do, but expect it to be SLOW. Better is to look for alternatives. For example, can you do at least one of those integrals analytically?

I would point out that far too often, I see people doing numerical integration of simple Gaussian density functions, when in fact, that is doable using a function call for the Gaussian CDF (or erfc if you wish to do the transformation yourself.)

Another choice is the use of Monte Carlo integration, which becomes relatively more efficient when you go into a higher number of dimensions.

[Mike Hosea]

True (curse of dimensionality). The integralN function is cranking... I have no idea when it will stop. :)

[John D'Errico]

Long ago, in galaxy far, far away, I did write a 5-fold integration code for a specific problem. Ran like a dog with two legs. The only thing I learned from that effort was to never try to solve that class of problem with brute force integration. We did not need the accuracy I thought I had to achieve, and there would have been many alternative routes I'd try now if I truly understood that fact. So, I suppose I did get something useful out of it, even though the project was a complete failure otherwise.

[sun]

Thanks for your answer here. John and Mike, you guys are very professional. When you are talking about

'sin(x) over the interval [0, 10*pi], took 390 function evaluations.'  and ' exp, over the interval [0, 2] took 150 function evaluations'

could you give me a link for me learn about this? Be honest, I don't know about dimensionality or function evaluations for matlab. Thank you again.

[Mike Hosea]

This is not really a "MATLAB" topic. John's pointing out that the amount of work required to do a 5-fold integral numerically via a quadrature rule is roughly the amount of work raised to the 5th power of doing a single integral. Sometimes that's a pessimistic analysis, but sometimes it's quite accurate. If you look at the documentation that I wrote in integralN, you'll see that I don't really recommend the approach. It's worth a try, but I think I used the words "excruciatingly slow" for 6-fold integrals. Your problem appears to be such in 5-fold integrals. I will see if it was successful tomorrow morning.

[sun]

Thank you very much, mike.I will check here back later.

[Mike Hosea]

The integration was unsuccessful:

>> myscript
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test. 
> In integral2Calc>funfun\private\integral2Calc>integral2t (line 129)
[snip]
q =
   NaN

Since the number of function evaluations was exceeded in a double integral sub-calculation, the boundary singularity may be a bit too strong for these codes to handle.

[John D'Errico]

I was a bit surprised, that the number of function evals taken was not apparently one of the things that would be returned by an integration tool, as one of the later arguments. (Any code that I would write would offer that as important and useful information, but I'm a bit over the top sometimes in terms of how complete I make my code.)

It is trivial to learn however. I simply wrote a kernel function for those problems that counted the number of times it was called. In fact, integral only called the kernel function once for exp(x), but it sent in 150 points in that call to be evaluated. the sin(x) kernel was a bit more curvy, it had two calls to the objective, with 390 total points to evaluate.

[sun]

Sad... Thanks alot, Mike and John. I am doing more research about it. I will be back later on.

[Mike Hosea]

I decided not to return the number of function evaluations because I think it is more confusing to the average user than it is helpful in the context of vectorized methods. Do we return 1 or 150? The former is literally correct, the latter is the number that you might connect with your past uses of this measure of "work". The 150 number, is, of course, just reflecting the conservative initial mesh that INTEGRAL uses. Possibly it could have gotten away with many fewer evaluations, and in fact, when INTEGRAL2 and INTEGRAL3 iterate with INTEGRAL's 1-D algorithm, they do use a coarser initial mesh.

What I have been wanting to do is develop some kind of visualization of the integration data that will help you to figure out what happened when integrations fail, something that generalizes well to all methods we might want to put under the "hood" of INTEGRAL, INTEGRAL2, and INTEGRAL3. My first stab at this was the "FailurePlot" in QUAD2D, but that was "free", and applying something like that for the 'iterated' method would require extra work and storage (as INTEGRAL does not maintain a master list of all abscissa's used), not to mention require some creative choices when limits are infinite. At the present time I recognize that there is a dearth of diagnostic information available.

Another thing I want is an "expert" interface to the integration routines that exposes all manner of knobs and switches as well as provides extra output data, such as the global error estimate and perhaps a lot more.