Problem with definite Integral

Question

AVM 2020-6-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/540251-problem-with-definite-integral

评论： AVM 2020-6-7

I am trying to solve the following definite double integration numerically. The expressions contain summaitions also but that is being executed within few seconds. When the double integration section comes, it is taking extremy huge time even after 7 hours it is still going on without any output. Any advice will be highly appreciated.

clc;
syms  n r theta m p  l t
w=1.0;
d=1.0;
g=0.2;
lmd=0.5;
assume(r,'real');
assume(theta,'real');
assume(t, 'real');
om=sqrt(((w).^2)-(4.*(g.^2)));
mu=sqrt((w+om)./(2.*om));
nu=((w-om)./(2.*g)).*mu;
eta=(((lmd)./((2.*g)+w)).*(1+((w-om)./(2.*g)))).*mu;
En=((n+(1./2)).*om)-(w./2)-(((lmd).^2)/((2.*g)+w));
um=((m+(1./2)).*om)-(w./2)-(((lmd).^2)/((2.*g)+w));
Enn=((n-(1./2)).*om)-(w./2)-(((lmd).^2)/((2.*g)+w));
umm=((m-(1./2)).*om)-(w./2)-(((lmd).^2)/((2.*g)+w));
Dn=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL(n,(4.*((eta).^2))));     
Dm=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL(m,(4.*((eta).^2))));  
Dnn=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL((n-1),(4.*(eta.^2)))); 
Dmm=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL((m-1),(4.*(eta.^2)))); 
Em=En - Dn;
Um=um-Dm;
Ep=Enn + Dnn;
Up=umm +Dmm;
epsn=(Ep-Em)./2; 
epsm=(Up-Um)./2;
Deln=(eta.*d./sqrt(n)).*exp(-2.*(eta.^2)).*laguerreL((n-1),1,(4.*(eta.^2)));
Delm=(eta.*d./sqrt(m)).*exp(-2.*(eta.^2)).*laguerreL((m-1),1,(4.*(eta.^2)));
xn=sqrt(((epsn).^2)+((Deln).^(2))); 
xm=sqrt(((epsm).^2)+((Delm).^(2))); 
zetapn=sqrt(((xn)+(epsn))./(2.*xn));
zetamn=sqrt(((xn)-(epsn))./(2.*xn));
zetapm=sqrt(((xm)+(epsm))./(2.*xm));
zetamm=sqrt(((xm)-(epsm))./(2.*xm));
z= 1i.*(mu-nu).*eta./(sqrt(2.*mu.*nu));
a1n=(zetapn./sqrt(factorial(n-1))).*((-nu./(2.*mu)).^(-1./2)).*hermiteH(n-1, z);
b1n=(Deln./abs(Deln)).*(zetamn./sqrt(factorial(n))).*hermiteH(n, z);
a2n=(zetamn./sqrt(factorial(n-1))).*((-nu./(2.*mu)).^(-1./2))*hermiteH(n-1, z);
b2n= (Deln./abs(Deln)).*(zetapn./sqrt(factorial(n))).*hermiteH(n, z);
a1m=(zetapm./sqrt(factorial(m-1))).*((-nu./(2.*mu)).^(-1./2)).*hermiteH(m-1, z);
b1m=(Delm./abs(Delm)).*(zetamm./sqrt(factorial(m))).*hermiteH(m, z);
a2m=(zetamm./sqrt(factorial(m-1))).*((-nu./(2.*mu)).^(-1./2))*hermiteH(m-1, z);
b2m= (Delm./abs(Delm)).*(zetapm./sqrt(factorial(m))).*hermiteH(m, z);
c0= -(1./sqrt(2.*mu)).*exp(-((eta.^2)./2)+ ((nu.*(eta).^2)./(2.*mu)));
cpn= -c0.*((-nu./(2.*mu)).^(n./2)).*(a1n - b1n);
cmn= -c0.*((-nu./(2.*mu)).^(n./2)).*(a2n + b2n);
cpm= -c0.*((-nu./(2.*mu)).^(m./2)).*(a1m - b1m);
cmm= -c0.*((-nu./(2.*mu)).^(m./2)).*(a2m + b2m);
E0=(om./2)-(w./2)-(((lmd).^2)./((2.*g)+w));
eg= E0-((d./2).*(exp(-2.*((eta).^(2)))));
ep=(1./2).*(Ep+ Em + (sqrt(((Ep-Em).^2)+(4.*((Deln).^2)))));
em=(1./2).*(Ep+ Em - (sqrt(((Ep-Em).^2)+(4.*((Deln).^2)))));
upp= (1./2).*(Up+ Um + (sqrt(((Up-Um).^2)+(4.*((Delm).^2)))));
umm= (1./2).*(Up+ Um - (sqrt(((Up-Um).^2)+(4.*((Delm).^2)))));
c0t= c0.*exp(-1i.*eg.*t);
cpnt= cpn.*exp(-1i.*ep.*t);
cmnt= cmn.*exp(-1i.*em.*t);
cpmt= cpm.*exp(-1i.*upp.*t);
cmmt= cmm.*exp(-1i.*umm.*t);
Ant=zetapn.*cpnt + zetamn.*cmnt;
Bnt= (Deln./abs(Deln)).*(zetamn.*cptn - zetapn.*cmnt);
Amt= zetapm.*cpmt + zetamm.*cmmt;
Bmt= (Delm./abs(Delm)).*(zetamm.*cpmt - zetapm.*cmmt);
beta= r.*exp(1i.*theta);
guard_digits = 10;
sp11= ((1i.^p)./factorial(p)).*((nu./(2.*mu)).^(p./2)).*hermiteH(p, 1i.*beta./sqrt(2.*mu.*nu)).*(eta.^(p+m)).*hypergeom([-p -m],[], -1./(eta.^2));
Hp11= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-((beta.^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(m))).*sum(vpa(subs(sp11,p,1:20), guard_digits));
sp22= (((-1i).^l)./factorial(l)).*((nu./(2.*mu)).^(l./2)).*hermiteH(l, -1i.*conj(beta)./sqrt(2.*mu.*nu)).*(eta.^(l+n)).*hypergeom([-l -n],[], -1./(eta.^2));
Hp22= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-(((conj(beta)).^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(n))).*sum(vpa(subs(sp22,l,1:20), guard_digits));
sm11= ((1i.^p)./factorial(p)).*((nu./(2.*mu)).^(p./2)).*hermiteH(p, 1i.*beta./sqrt(2.*mu.*nu)).*(-eta.^(p+m)).*hypergeom([-p -m],[], -1./(eta.^2));
Hm11= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-((beta.^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(m))).*sum(vpa(subs(sm11,p,1:20), guard_digits));
sm22= (((-1i).^l)./factorial(l)).*((nu./(2.*mu)).^(l./2)).*hermiteH(l, -1i.*conj(beta)./sqrt(2.*mu.*nu)).*(-eta.^(l+n)).*hypergeom([-l -n],[], -1./(eta.^2));
Hm22= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-(((conj(beta)).^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(n))).*sum(vpa(subs(sm22,l,1:20), guard_digits));
Hp1=Hp22.*Hp11;
Hm1=Hm22.*Hm11;
Hp(n,m)= (1./(2.*pi)).*(Hp1 + Hm1);
Hm(n,m)= (1./(2.*pi)).*(Hp1 - Hm1);
f11=((abs(c0t)).^2).*Hp(0,0);
f22= c0t.*conj(Ant).*Hm(0,n-1) + conj(c0t).*Ant.*Hm(n-1,0)+ c0t.*conj(Bnt).*Hp(0,n) + conj(c0t).*Bnt.*Hp(n,0);
f33= Ant.*conj(Amt).*Hp(n-1,m-1) + Bnt.*conj(Bmt).*Hp(n,m) + Bnt.*conj(Amt).*Hm(n,m-1) +Ant.*conj(Bmt).*Hm(n-1,m);
sf33= sum(vpa(subs(f33,m,1:20), guard_digits));
f=f11 + sum(vpa(subs(f22,n,1:20), guard_digits)) + sum(vpa(subs(sf33,n,1:20), guard_digits));
vpaintegral(vpaintegral(f, r, [0 10]), theta, [0 2.*pi]) %% 'r' and 'theta' are integration variable
%int(int(f,r,0,10),theta,0,2*pi) 

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Ameer Hamza 2020-6-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/540251-problem-with-definite-integral#answer_444847

在 MATLAB Online 中打开

If you check at the expression of 'f', you can see it also has 't'. So even if you try to numerically integrate it w.r.t. 'r' and 'theta', the answer will still be symbolic.

Also, I suggest you to use matlabFunction() to convert the symbolic expression into a floating-point function, which is much faster than the symbolic calculations. For example, instead of vpaintegral(), try this

F = matlabFunction(f, 'Vars', [r theta t]);
int_val = integral2(@(r, theta) F(r, theta, 0), 0, 10, 0, 2*pi)

This assumes that t=0 to get a function in terms of r and theta.

15 个评论
显示 13更早的评论隐藏 13更早的评论

AVM 2020-6-2

编辑：AVM 2020-6-7

在 MATLAB Online 中打开

@Ameer: After rigorous checking of my equations , the current version of it given below. Now the problem is Long-run time even though this is the correct equations. Here, I put t=0 ((Actually I have to chose some 't' according to my system of Eqns. )) but still it is taking so much time to execute as the probelm was in begin with. Although upto the function ''f'' it is executing within a min. Please, have a look and advise me what should I do.

clc;
syms  n r theta m p  l 
w=1.0;
d=0.3;
g=0.0001;
lmb=0.1;
t=0;
assume(r,'real');
assume(theta,'real');
om=sqrt(((w).^2)-(4.*(g.^2)));
mu=sqrt((w+om)./(2.*om));
nu=((w-om)./(2.*g)).*mu;
eta=(((lmb)./((2.*g)+w)).*(1+((w-om)./(2.*g)))).*mu;
En=((n+(1./2)).*om)-(w./2)-(((lmb).^2)/((2.*g)+w));
Enn=((n-(1./2)).*om)-(w./2)-(((lmb).^2)/((2.*g)+w));
engpm=((m+(1./2)).*om)-(w./2)-(((lmb).^2)/((2.*g)+w));
engmm=((m-(1./2)).*om)-(w./2)-(((lmb).^2)/((2.*g)+w));
Dn=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL(n,(4.*((eta).^2))));    
Dm=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL(m,(4.*((eta).^2))));  
Dnn=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL((n-1),(4.*(eta.^2)))); 
Dmm=(d./2).*(exp(-2.*((eta).^(2)))).*(laguerreL((m-1),(4.*(eta.^2)))); 
Em=En - Dn;
Ep=Enn + Dnn;
Um=engpm-Dm;
Up=engmm +Dmm;
epsn=(Ep-Em)./2; 
epsm=(Up-Um)./2;
Deln=(eta.*d./sqrt(n)).*exp(-2.*(eta.^2)).*laguerreL((n-1),1,(4.*(eta.^2)));
Delm=(eta.*d./sqrt(m)).*exp(-2.*(eta.^2)).*laguerreL((m-1),1,(4.*(eta.^2)));
xn=sqrt(((epsn).^2)+((Deln).^(2)));
xm=sqrt(((epsm).^2)+((Delm).^(2))); 
zetapn=sqrt(((xn)+(epsn))./(2.*xn));
zetamn=sqrt(((xn)-(epsn))./(2.*xn));
zetapm=sqrt(((xm)+(epsm))./(2.*xm));
zetamm=sqrt(((xm)-(epsm))./(2.*xm));
z= 1i.*(mu-nu).*eta./(sqrt(2.*mu.*nu));
a1n=(zetapn./sqrt(factorial(n-1))).*((-nu./(2.*mu)).^(-1./2)).*hermiteH(n-1, z);
b1n=(Deln./abs(Deln)).*(zetamn./sqrt(factorial(n))).*hermiteH(n, z);
a2n=(zetamn./sqrt(factorial(n-1))).*((-nu./(2.*mu)).^(-1./2))*hermiteH(n-1, z);
b2n= (Deln./abs(Deln)).*(zetapn./sqrt(factorial(n))).*hermiteH(n, z);
a1m=(zetapm./sqrt(factorial(m-1))).*((-nu./(2.*mu)).^(-1./2)).*hermiteH(m-1, z);
b1m=(Delm./abs(Delm)).*(zetamm./sqrt(factorial(m))).*hermiteH(m, z);
a2m=(zetamm./sqrt(factorial(m-1))).*((-nu./(2.*mu)).^(-1./2))*hermiteH(m-1, z);
b2m= (Delm./abs(Delm)).*(zetapm./sqrt(factorial(m))).*hermiteH(m, z);
c0= -(1./sqrt(2.*mu)).*exp(-((eta.^2)./2)+ ((nu.*(eta).^2)./(2.*mu)));
cpn= -c0.*((-nu./(2.*mu)).^(n./2)).*(a1n - b1n);
cmn= -c0.*((-nu./(2.*mu)).^(n./2)).*(a2n + b2n);
cpm= -c0.*((-nu./(2.*mu)).^(m./2)).*(a1m - b1m);
cmm= -c0.*((-nu./(2.*mu)).^(m./2)).*(a2m + b2m);
guard_digits = 10;
E0=(om./2)-(w./2)-(((lmb).^2)./((2.*g)+w));
eg= E0-((d./2).*(exp(-2.*((eta).^(2)))));
ep=(1./2).*(Ep+ Em + (sqrt(((Ep-Em).^2)+(4.*((Deln).^2)))));
em=(1./2).*(Ep+ Em - (sqrt(((Ep-Em).^2)+(4.*((Deln).^2)))));
upp= (1./2).*(Up+ Um + (sqrt(((Up-Um).^2)+(4.*((Delm).^2)))));
umm= (1./2).*(Up+ Um - (sqrt(((Up-Um).^2)+(4.*((Delm).^2)))));
c0t= c0.*exp(-1i.*eg.*t);
cpnt= cpn.*exp(-1i.*ep.*t);
cmnt= cmn.*exp(-1i.*em.*t);
cpmt= cpm.*exp(-1i.*upp.*t);
cmmt= cmm.*exp(-1i.*umm.*t);
Ant=zetapn.*cpnt + zetamn.*cmnt;
Bnt= (Deln./abs(Deln)).*(zetamn.*cpnt - zetapn.*cmnt);
Amt= zetapm.*cpmt + zetamm.*cmmt;
Bmt= (Delm./abs(Delm)).*(zetamm.*cpmt - zetapm.*cmmt);
beta= r.*exp(1i.*theta);
sp11= ((1i.^p)./factorial(p)).*((nu./(2.*mu)).^(p./2)).*hermiteH(p, 1i.*beta./sqrt(2.*mu.*nu)).*(eta.^(p+m)).*hypergeom([-p -m],[], -1./(eta.^2));
Hp11= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-((beta.^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(m))).*sum(vpa(subs(sp11,p,1:5), guard_digits));
sp22= (((-1i).^l)./factorial(l)).*((nu./(2.*mu)).^(l./2)).*hermiteH(l, -1i.*conj(beta)./sqrt(2.*mu.*nu)).*(eta.^(l+n)).*hypergeom([-l -n],[], -1./(eta.^2));
Hp22= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-(((conj(beta)).^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(n))).*sum(vpa(subs(sp22,l,1:5), guard_digits));
sm11= ((1i.^p)./factorial(p)).*((nu./(2.*mu)).^(p./2)).*hermiteH(p, 1i.*beta./sqrt(2.*mu.*nu)).*((-eta).^(p+m)).*hypergeom([-p -m],[], -1./(eta.^2));
Hm11= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-((beta.^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(m))).*sum(vpa(subs(sm11,p,1:5), guard_digits));
sm22= (((-1i).^l)./factorial(l)).*((nu./(2.*mu)).^(l./2)).*hermiteH(l, -1i.*conj(beta)./sqrt(2.*mu.*nu)).*((-eta).^(l+n)).*hypergeom([-l -n],[], -1./(eta.^2));
Hm22= ((exp(-((eta.^2)./2)-(((abs(beta)).^2)./2)-(((conj(beta)).^2).*(nu)./(2.*mu))))./sqrt(mu.*factorial(n))).*sum(vpa(subs(sm22,l,1:5), guard_digits));
Hp1=Hp22.*Hp11;
Hm1=Hm22.*Hm11;
Hp(n,m)= (1./(2.*pi)).*(Hp1 + Hm1);
Hm(n,m)= (1./(2.*pi)).*(Hp1 - Hm1);
f11=((abs(c0t)).^2).*Hp(0,0);
f22= c0t.*conj(Ant).*Hm(0,n-1) + conj(c0t).*Ant.*Hm(n-1,0)+ c0t.*conj(Bnt).*Hp(0,n) + conj(c0t).*Bnt.*Hp(n,0);
f33= Ant.*conj(Amt).*Hp(n-1,m-1) + Bnt.*conj(Bmt).*Hp(n,m) + Bnt.*conj(Amt).*Hm(n,m-1) + Ant.*conj(Bmt).*Hm(n-1,m);
sf33= sum(vpa(subs(Q33,m,1:5), guard_digits));
f= f11 + sum(vpa(subs(f22,n,1:5), guard_digits)) + sum(vpa(subs(sf33,n,1:5), guard_digits));
F = matlabFunction(f, 'Vars', [r theta]);
int_val = integral2(@(r, theta) F(r, theta), 0, 5, 0, 2*pi)
polarplot(theta, int_val,'r')

AVM 2020-6-7

Thanks for your advice. By the way I have reduced the number of summations involved in the equations ( upto 1-5) and now it is running much faster than previous. Now I need to plot the ''int_val'' w.r.t. '' theta'' in polar plot without using ''int(int())'' command because when I need to plot that w.r.t. ''theta'', I think ''Integral2''will no longer work. So what should I do in this case to get the plot.

I have aready used ''int(int())'' but it is taking extremely long time. Pl help me to solve this problem

AVM 2020-6-7