MATLAB: Problem with definite Integral

double integration

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)

Best Answer

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.
Related Question