speed up infinite integral with 4 variables

Hi;

I have an infinite integral with 4 variables and i coded as follows. But it takes too much time. Matlab is working but there is no answer. when i changed abstol and reltol integration yields unexpected results. Is there any way to speed up this integration with high tolerances?

My codes are as;

function Fig1_pc_Gauss_scint_vs_norm_source_size 
clc

global s1x;
global s1y;
global s2x;
global s2y;

N=4;    %beam number
alphas=0.05;    %beam source size
L=4000;    %distance
l0=0.001;  %inner scale 1 mm
L0=25;    %outer scale 25 m
Cn2=1e-15;      %turbulence structure constant
lambda=1.55e-6;  %wavelength
k=2*pi/lambda;

ksi1=1;    %anisotropic factor
ksi2=3;
ksi3=5;
ksi4=8;
ksi5=10;

a=-0.061;
b=2.836;

kpp0=2*pi/L0;

j=1;
I0L1=0;I0L2=0;I0L3=0;I0L4=0;I0L5=0;
I0V=0;%I0V2=0;I0V3=0;I0V4=0;I0V5=0;I0V6=0; 

%inl = alphas*100;
xmin = -inf; xmax = inf; ymin = -inf; ymax = inf; zmin = -inf; zmax = inf; wmin = -inf; wmax = inf;

for alpha=3.1:0.1:3.9

     alphan(j)=alpha;

     Aalpha=gamma(alpha-2)*sin(pi*(alpha-3)/2)/(4*pi*pi);
     c1alpha=(pi*Aalpha*(gamma(3/2-alpha/2)*((3-alpha)/3)+a*gamma(2-alpha/2)*((4-
     alpha)/3)+b*gamma(3-(3*alpha)/4)*((4-alpha)/2)))^(1/(alpha-5));
     kppH=c1alpha/l0;

     n_alfa = 0.5*gamma(alpha)*((2*pi)^(-11/6+(alpha)/2))*(((lambda*(L))^(11/6-(alpha)/2)));
     d_alfa = gamma(1-alpha/2)*((gamma(alpha/2))^2)*gamma(alpha-1)*cos(pi*alpha/2)*sin(pi*alpha/4);
     Cn2tilda = (-1)*(n_alfa./d_alfa)*Cn2;

       ro01=(1.629*k*k*L*Aalpha*Cn2tilda*(ksi1^(-2))*((kpp0^(4-alpha))*hypergeom(2,7/6,(kpp0*kpp0/(kppH*kppH)))+...
       a*(kpp0^(5-alpha))*hypergeom(5/2,5/3,(kpp0*kpp0/(kppH*kppH)))/kppH+...
       b*(kpp0^(7-(3*alpha)/2))*gamma((14-alpha)/4)*hypergeom((14-alpha)/4,(32-3*alpha)/12,
       (kpp0*kpp0/(kppH*kppH)))/(kppH^(3-alpha/2))))^(-0.5);

       ro02=(1.629*k*k*L*Aalpha*Cn2tilda*(ksi2^(-2))*((kpp0^(4-alpha))*hypergeom(2,7/6,(kpp0*kpp0/(kppH*kppH)))+...
       a*(kpp0^(5-alpha))*hypergeom(5/2,5/3,(kpp0*kpp0/(kppH*kppH)))/kppH+...
       b*(kpp0^(7-(3*alpha)/2))*gamma((14-alpha)/4)*hypergeom((14-alpha)/4,(32-3*alpha)/12,
       (kpp0*kpp0/(kppH*kppH)))/(kppH^(3-alpha/2))))^(-0.5);

       ro03=(1.629*k*k*L*Aalpha*Cn2tilda*(ksi3^(-2))*((kpp0^(4-alpha))*hypergeom(2,7/6,(kpp0*kpp0/(kppH*kppH)))+...
       a*(kpp0^(5-alpha))*hypergeom(5/2,5/3,(kpp0*kpp0/(kppH*kppH)))/kppH+...
       b*(kpp0^(7-(3*alpha)/2))*gamma((14-alpha)/4)*hypergeom((14-alpha)/4,(32-3*alpha)/12,
       (kpp0*kpp0/(kppH*kppH)))/(kppH^(3-alpha/2))))^(-0.5);

       ro04=(1.629*k*k*L*Aalpha*Cn2tilda*(ksi4^(-2))*((kpp0^(4-alpha))*hypergeom(2,7/6,(kpp0*kpp0/(kppH*kppH)))+...
       a*(kpp0^(5-alpha))*hypergeom(5/2,5/3,(kpp0*kpp0/(kppH*kppH)))/kppH+...
       b*(kpp0^(7-(3*alpha)/2))*gamma((14-alpha)/4)*hypergeom((14-alpha)/4,(32-3*alpha)/12,
       (kpp0*kpp0/(kppH*kppH)))/(kppH^(3-alpha/2))))^(-0.5);

       ro05=(1.629*k*k*L*Aalpha*Cn2tilda*(ksi5^(-2))*((kpp0^(4-alpha))*hypergeom(2,7/6,(kpp0*kpp0/(kppH*kppH)))+...
       a*(kpp0^(5-alpha))*hypergeom(5/2,5/3,(kpp0*kpp0/(kppH*kppH)))/kppH+...
       b*(kpp0^(7-(3*alpha)/2))*gamma((14-alpha)/4)*hypergeom((14-alpha)/4,(32-3*alpha)/12,
       (kpp0*kpp0/(kppH*kppH)))/(kppH^(3-alpha/2))))^(-0.5)

for l1=1:1:N

      alphasl1=alphas/sqrt(l1);
      alphal1=1/(2*k*alphasl1);

      Al1=(((-1)^(l1-1))*factorial(N)/(factorial(l1)*factorial(N-l1)*N));
    
      for l2=1:1:N
    
          alphasl2=alphas/sqrt(l2);
          alphal2=1/(2*k*alphasl2);
    Al2=(((-1)^(l2-1))*factorial(N)/(factorial(l2)*factorial(N-l2)*N));
    
    fm1 =@(s1x,s1y,s2x,s2y)(exp((-k*alphal1+1i*k/(2*L))*s1x.*s1x).*exp((-k*alphal1+1i*k/(2*L))*s1y.*s1y).*exp((-k*alphal2-1i*k/(2*L))*s2x.*s2x).*...
         exp((-k*alphal2-1i*k/(2*L))*s2y.*s2y).*((((s1x-s2x).*(s1x-s2x)+(s1y-s2y).*(s1y-s2y)).^((alpha-2)/2))*(-1)/(ro01^(alpha-2))));
    tcm1 = integralN(fm1,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax)
    
    fm2 =@(s1x,s1y,s2x,s2y)(exp((-k*alphal1+1i*k/(2*L))*s1x.*s1x).*exp((-k*alphal1+1i*k/(2*L))*s1y.*s1y).*exp((-k*alphal2-1i*k/(2*L))*s2x.*s2x).*...
        exp((-k*alphal2-1i*k/(2*L))*s2y.*s2y).*((((s1x-s2x).*(s1x-s2x)+(s1y-s2y).*(s1y-s2y)).^((alpha-2)/2))*(-1)/(ro02^(alpha-2))));
    tcm2 = integralN(fm2,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax)
    
    fm3 =@(s1x,s1y,s2x,s2y)(exp((-k*alphal1+1i*k/(2*L))*s1x.*s1x).*exp((-k*alphal1+1i*k/(2*L))*s1y.*s1y).*exp((-k*alphal2-1i*k/(2*L))*s2x.*s2x).*...
    exp((-k*alphal2-1i*k/(2*L))*s2y.*s2y).*((((s1x-s2x).*(s1x-s2x)+(s1y-s2y).*(s1y-s2y)).^((alpha-2)/2))*(-1)/(ro03^(alpha-2))));
    tcm3 = integralN(fm3,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax)
    
    fm4 =@(s1x,s1y,s2x,s2y)(exp((-k*alphal1+1i*k/(2*L))*s1x.*s1x).*exp((-k*alphal1+1i*k/(2*L))*s1y.*s1y).*exp((-k*alphal2-1i*k/(2*L))*s2x.*s2x).*...
    exp((-k*alphal2-1i*k/(2*L))*s2y.*s2y).*((((s1x-s2x).*(s1x-s2x)+(s1y-s2y).*(s1y-s2y)).^((alpha-2)/2))*(-1)/(ro04^(alpha-2))));
    tcm4 =  integralN(fm4,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax)
    
    fm5 =@(s1x,s1y,s2x,s2y)(exp((-k*alphal1+1i*k/(2*L))*s1x.*s1x).*exp((-k*alphal1+1i*k/(2*L))*s1y.*s1y).*exp((-k*alphal2-1i*k/(2*L))*s2x.*s2x).*...
    exp((-k*alphal2-1i*k/(2*L))*s2y.*s2y).*((((s1x-s2x).*(s1x-s2x)+(s1y-s2y).*(s1y-s2y)).^((alpha-2)/2))*(-1)/(ro05^(alpha-2))));
          tcm5 = integralN(fm5,xmin,xmax,ymin,ymax,zmin,zmax,wmin,wmax)
    
    I0L1=I0L1+Al1*Al2*tcm1;
          I0L2=I0L2+Al1*Al2*tcm2;
       I0L3=I0L3+Al1*Al2*tcm3;
    I0L4=I0L4+Al1*Al2*tcm4;
    I0L5=I0L5+Al1*Al2*tcm5;
       

  I0V=I0V+Al1*Al2/(((k*alphal1-(1i)*k/(2*L)))*(k*alphal2+(1i)*k/(2*L)));
  
    end

end

trns1(j)=I0L1/I0V

trns2(j)=I0L2/I0V

trns3(j)=I0L3/I0V

trns4(j)=I0L4/I0V

trns5(j)=I0L5/I0V

j=j+1;

end

hold on

plot(alphan,trns1,'--k','LineWidth',2);

plot(alphan,trns2,'-k','LineWidth',2);

plot(alphan,trns3,':k','LineWidth',2);

plot(alphan,trns4,'-.k','LineWidth',2);

plot(alphan,trns5,'--b','LineWidth',2);

hold off

xlabel( '\it\alpha\rm','FontSize',24,'FontSize',24,'FontName','Times New Roman'); set(gca,'FontSize',24,'FontName','Times New Roman');

ylabel('Transmittance, \tau','FontSize',24,'FontName','Times New Roman');set(gca,'FontSize',24,'FontName','Times New Roman');

[Walter Roberson]

I find that I have difficulty understanding the reasoning for wanting integration to "high tolerance" when you cannot decide how many significant figures your coefficients are. For example 0.05 is one significant figure, but 2.836 is four significant figures. With only one significant figure for your coefficients, why would you be concerned with anything other than the "order of magnitude" of the integral ?

And remember,

for alpha=3.1:0.1:3.9

is not the same thing numerically as

for alpha_t=31:1:39
  alpha = alpha_t / 10;

> mat2str((3.1:0.1:3.9) - (31:1:39)/10)
[0 0 4.44089209850063e-16 4.44089209850063e-16 0 -4.44089209850063e-16 -4.44089209850063e-16 0 0]

[ylcnt]

do you mean ı change the tolerances and limits?

[Walter Roberson]

"Ray's Rule of Precision: Measure with a micrometer. Mark with chalk. Cut with an axe."

You are more trying to measure with an axe and cut with a micrometer. Your coefficients at the moment cannot justify asking for high precision output, because you do not have high precision input.

What is it that you are calculating? Some kind of multimode fibre transmission?

[ylcnt]

This is an optical turbulence calculation. There is no analytical solution.

Answer 1

[John D'Errico]

Sorry, but you do know what ROFLMAO means?

Big problems take big time. Computers are not infinitely fast.

Check the time required for any single function evaluation. A 4-fold numerical integration will certainly require probably millions of function evals, IF you are lucky. I'd be not at all surprised if you should be expecting something at least on the order of 1e8 or 1e9 function evals here. And you want high accuracy, so maybe more. Look at it like this, a 1-d application of an adaptive numerical integration will require on the order of 100-200 kernel evals, on any non-trivial function. 100^4 = 1e8.

Again, ROFLMAO. Yeah, I know, nobody likes having their problems being laughed at. But you need to understand the issues.

You might consider limiting the domain of integration. Infinity is a long way off, and if your kernel is effectively zero, then why bother looking out there in the hinterlands?

Were this my problem to solve, I'd start by looking at the form of your kernel, and decide if some variety of Gaussian integration might apply. Your code is too messy for me to even bother with that. Of course, you would need to choose an appropriate weight function, so the proper family of Gaussian quadrature would be important.