Problem with definite Integral

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) 

Answer 1

[Ameer Hamza]

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.

[AVM]

@Ameer: Thanks for your advice. Now, I am trying as you said.

[AVM]

It is very nice. Now it is giving output within few mins. One thing that I would like to know when it assumes t=0, does it mean that t =0 through out the calculation? Pl let me know.

[Ameer Hamza]

Yes, it put t=0 throughout the calculation. Whether 't' should be present in the final expression depends on the system of the equation you are trying to implement.

[AVM]

okay. So, even it assumes t =0 initially, it does not mean that 't' will be missing from the output expression ? Am I right? So, it depends upon the function that I am trying to integrate over. It may happen that in some cases 't' may present ultimately in the output depending upon equations. Am I correct? Pl help me to understand this point.

[Ameer Hamza]

Do you want 't' to be present in the final expression? Is it the correct way according to your system of equations?

[AVM]

Actually accroding to my system of equations, the 't' will not present ultimately in the output after integration. In the process of entire evolution 't' will be somehow vanished. That is my goal that output will be independent of 't'. But I am asking this question that is this happening (vanish of 't' from output) due to the assumption that 't =0' initially? I would like to understand this point only. Pl make me clear.

[AVM]

Here, also I am getting the output independent of 't' which I am expecting. But Is this happening only for the underlying assumption that 't =0' ?

[Ameer Hamza]

My guess is that there is some mistake in writing the system of equation. I also think 't' should not be present in the expression of 'f'. You may need to check if all equations are written correctly and find out why 't' is not vanishing in the final expression of 'f'.

The way I am doing it (putting t=0) does not seems to be a correct solution to me. I just gave it as an example to show how integral2 works. From your question, it appears that the expression of 'f' should only contain r and theta.

Note that t=0 is just an arbitrary choice. You can choose t=1 (or any other constant value), and the integral will still work.

[AVM]

Thanks for your clarification. I am now checking the expressions.

[AVM]

@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')

[Ameer Hamza]

I am not sure why it is taking such a long time now. Maybe keep 't' as variable as you did in previous code and then put t=0 after conversion to function handle as I did in my answer.

[AVM]

I have tried in both ways putting 't' =0 initially as well as t =0 later on as you have shown in your code. Unfortunately in both cases it is taking long time. Is there any other alternative command to resolve the issue?

[Ameer Hamza]

This is caused by the complexity of your system of equations. I don't think there is a way to improve the calculation further in MATLAB. It will need to decrease the complexity of your equations somehow to get any speed improvement.

[AVM]

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]

@Ameer: Okay the issue is not relevant here. I just asked this question by mistake. It will no longer releted to my case now.