How to do numerical 4-fold integral (more than triple) by using matlab?

Dear friend, My question is that how to do the 4-fold integration. I use to want to do more that 5-fold, but I failed..... so, I decide to go back to a lower level, just 4 fold numerical integration. Please check the code below, there are some parameters and one function.

clear all; 
%%%=== 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
  powe= 0 ; 
  seta= 10^powe;
  q=pi.*(la1.*p1.^(2./a)+la2.*p2.^(2./a));
%%%=== just some parameters END here ==============  

%%%=== Function I want to do 4-fold numerical integration below ===

    yi =   exp(-u3.*(1+2.*...
        (  out1      )./...
    (  (( (u3./u1).^(a./2) + (u3./u2).^(a./2) + 1    ).^(2./a)).*seta.^(-2./a)))).*...
                  exp(-sgma2.*q.^(-a./2).* seta.*u3.^(a./2)./...
          ((( (u3./u1).^(a./2) + (u3./u2).^(a./2) + 1    ).^(2./a)).^(a./2))  );

Ok, As you seen, variables are u1,u2,u3,u4. There will be 4 'nest' integration. The most inside variable is u1, 2nd is u2, 3rd is u3, the most outside is u4. The range of u1 is [0, u2], u2 is [0, u3], u3 is [0, u4], u4=[0, Inf]. This is a numerical integration, this means after 4 level numerical integration, the final output is a real number, not a symbol.

Thank you so much for your time. Thanks. %%%%%====================================================================== The information below is what I've tried.

In another "try-version" of my code, yi is only function of u1, u2. and if I use int 2 times as below,

    %%= 1000 means Inf here, 1st int var is u1 and range is[0 u2], 2nd int var is w.r.t u2 and [0, 1000 or Inf]

out2 = double ( int(int(yi,u1,0,u2),u2, 0 ,1000))

For the code above, it's working! it gave me a real number. Then I want to try 3 fold, If I do the same thing for 3 fold integration, I failed. So, I change yi to function of u1,u2,u3. and I do code below,

    out2 = double ( int(int(int(yi,u1,0,u2),u2, 0 ,u3)),u3,0,inf)

1st var is u1, range is [0 u2]; 2nd var is u2, the range is [0 u3]; last one is u3, [0 Inf]. error from above code is

Warning: Explicit integral could not be found. 
Warning: Explicit integral could not be found. 
Warning: Explicit integral could not be found. 
Error using sym/double
Too many input arguments.

Error in ref7_equ11n2 (line 117)
out2 = double ( int(int(int(yi,u1,0,u2),u2, 0 ,u3)),u3,0,inf)

As you seen, I am not using quad2d, integral3 or other in-built functions. Because I know I will get a 'right' answer for 2 or 3 fold integration. But my target is for 4-fold. that's why I need to do it in 'int' way, step by step. Thank you.

[John D'Errico]

You seem to think that before you define a function of u1, u2, u3, u4, that you need to define them as syms. There is no need to do so. Just a waste of CPU cycles.

[John D'Errico]

Just a few off the cuff thoughts here late at night...

It looks like your region is a simplex in 4-dimensions, but with one leg that extends over [0,inf). So I would either think of doing the u4 integration in a Gaussian form (maybe Laguerre?), or I would choose some logical transformation from [0,inf) into some finite interval.

I might also think about refining that simplex into smaller pieces, then you could use Gaussian quadrature over those simplexes. It is late though for me, and my mind is toast.

[sun]

Dear John, could you say it more clear? I am not very follow you. So, which function is right one to do this 4-fold intergration? Thank you,

[sun]

Well, If I put a vpa before double, I will have the right answer for 3-fold integration. It's also working for 4-fold integration. thanks guys

Answer 1

[sun]

Well, If I put a vpa before double, I will have the right answer for 3-fold integration. It's also working for 4-fold integration. thanks guys