How to do numerical multiple integral (more than triple) by using matlab?
51 views (last 30 days)
Show older comments
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.
1 Comment
Mike Hosea
on 29 Dec 2014
You have failed to supply any arguments to out1fcn. Remember that function handles always require arguments when you are evaluating them.
Accepted Answer
Mike Hosea
on 28 Dec 2014
Edited: Mike Hosea
on 28 Dec 2014
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.
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)
3 Comments
Mike Hosea
on 29 Dec 2014
Edited: Mike Hosea
on 29 Dec 2014
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
on 30 Dec 2017
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.
More Answers (2)
Shoaibur Rahman
on 28 Dec 2014
Edited: Shoaibur Rahman
on 28 Dec 2014
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.
2 Comments
Shoaibur Rahman
on 29 Dec 2014
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)
John D'Errico
on 28 Dec 2014
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.
9 Comments
Mike Hosea
on 29 Dec 2014
Edited: Mike Hosea
on 29 Dec 2014
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.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!