A 4D Nested numerical integration

Hello. I am new to matlab. I am trying to integrate a function over time and 3-space coordinate using vpaintegral command given below by the S12int1 command. I have waited tens of hours but did not get any result. When it is over 2D given by S12int2 I get the answer in few seconds. In 3D I get a result in about 10 minutes. I list my questions below. Any help is very much appreciated.

1. What is the stragegy to evaluate multidimensional integrals in general ? Do you recommend vpaintegral over integrate command
2. In vpaintegral command I can set the accuracy by setting for instance 'RelTol',1e-15,'AbsTol',1e-20. What I observed is if the integrand's magnitude is much larger than the value set by the tolerance values above, integration takes forever. Why is that so ? Do I need to adjust tolerance values in accordance with the integrand's magnitude ?
3. When using nested vpaintegral commands Do I need to specify error tolerances for each vpaintegral command ? Or just for the inner most integral ? A related question is do I need to put the command 'ArrayValued',true in each nest ?
4. Is there any other approach ? I am condering turning integrals into summations and perform them in nested for loops. Do you think this would increase the error ?

Thank you!

close all;
%clc;
clear all;
syms t
syms phi
syms z
syms r


tic

zr=1/100; omega=1/100; tau=1; c=1; E0=1/10; kg=2; phik=pi/4; thetak=pi/7;

S12(t,r,z,phi) = exp(-1i*c*kg*t+1i*kg*z*cos(thetak)+ 1i*kg*r*cos(phi-phik)*sin(thetak)) ...
*exp(-(t+z/c)^2/tau^2)*exp(-2*(t-z/c)^2/tau^2)*zr^3/(z^2+zr^2)^(3/2)*exp(-3*r^2*zr*omega/(2*c*(z^2+zr^2)))...
*cos((z/c-t)*omega+r^2*z*omega/(2*c*(z^2+zr^2))-atan(z/zr))*cos((z/c-t)*omega+r^2*z*omega/(2*c*(z^2+zr^2))-atan(z/zr)) ...
*cos((-z/c-t)*omega-r^2*z*omega/(2*c*(z^2+zr^2))+atan(z/zr))*r;




%S12int1= vpaintegral(vpaintegral(vpaintegral(vpaintegral(S12(t,r,z,phi),t,[-3,3]),z,[-10,10]),r,[0,5]),phi,[0,2*pi]);

S12int2= vpaintegral(vpaintegral(vpaintegral(S12(t,r,z,pi/3),t,[-3,3]),r,[0,10]),z,[-10,10],'ArrayValued',true);


S12int2
%double(S12(1,1,1,1))




timeElapsed = toc

[Torsten]

The strategy is not to use a symbolic setup, but a pure numerical approach with integral/integral2/integral3. Further, as a first shot, use the default tolerances.

[Taqu3]

Hi Torsten, thank you for your comment. Is there integral4 version ? or is possible to use integral3 combined with integral ?

[Torsten]

Here is a person who already prepared this combination:

https://de.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

[Taqu3]

Is there a tutorial on how to use the file you have linked above ? Btw, I have tried the following

close all;

%clc;

clear all;

zr=1/100; omega=1/100; tau=1; c=1; E0=1/10; kg=2; phik=pi/4; thetak=pi/7;

S12 = @(t,r,z,phi) exp(-1i.*c.*kg.*t+1i.*kg.*z.*cos(thetak)+ 1i.*kg.*r.*cos(phi-phik).*sin(thetak)) ...

.*exp(-(t+z/c).^2/tau.^2).*exp(-2.*(t-z/c).^2/tau.^2).*zr.^3/(z.^2+zr.^2).^(3/2).*exp(-3.*r.^2.*zr.*omega/(2.*c.*(z.^2+zr.^2)))...

.*cos((z/c-t).*omega+r.^2.*z.*omega/(2.*c.*(z.^2+zr.^2))-atan(z/zr)).*cos((z/c-t).*omega+r.^2.*z.*omega/(2.*c.*(z^.2+zr.^2))-atan(z/zr)) ...

.*cos((-z/c-t).*omega-r.^2.*z.*omega/(2.*c.*(z.^2+zr.^2))+atan(z/zr)).*r;

%integral(@(phi)arrayfun(@(phi)integral3(@(t,r,z)S12(t,r,z,phi),-10,10,0,5,-10,10),phi),0.01,3,'AbsTol',1e-20);

integral(@(phi)integral3(@(t,r,z)S12(t,r,z,phi),-10,10,0,5,-10,10),0,3,'ArrayValued',true,'AbsTol',1e-20)

Both of them gave the warning messages such as

Warning: Matrix is singular to working precision.

Arrays have incompatible sizes for this operation.

[Torsten]

If you open integralN, you'll see some comments and examples from line 1 to line 75.

Start with simpler examples than yours to see how the code works.