speed up infinite integral with 4 variables

3 次查看（过去 30 天）
ylcnt
ylcnt 2017-5-12
评论： ylcnt 2017-5-12
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');
  4 个评论
显示 2更早的评论
Walter Roberson
Walter Roberson 2017-5-12
"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
ylcnt 2017-5-12
This is an optical turbulence calculation. There is no analytical solution.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

John D'Errico
John D'Errico 2017-5-12
编辑：John D'Errico 2017-5-12
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.

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by