Problem with definite Integral

AVM
2020 6 2
1 个回答
回答已采纳
更新时间：2020 6 7
3 次查看（30 天）

0 个投票

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)

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

 采纳的回答

Ameer Hamza
Ameer Hamza 2020-6-2

0 个投票

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
AVM 2020-6-2
@Ameer: Thanks for your advice. Now, I am trying as you said.
AVM
AVM 2020-6-2
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
Ameer Hamza 2020-6-2
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
AVM 2020-6-2
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
Ameer Hamza 2020-6-2
Do you want 't' to be present in the final expression? Is it the correct way according to your system of equations?
AVM
AVM 2020-6-2
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
AVM 2020-6-2
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
Ameer Hamza 2020-6-2
编辑：Ameer Hamza 2020-6-2
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
AVM 2020-6-2
Thanks for your clarification. I am now checking the expressions.
AVM
AVM 2020-6-2
编辑：AVM 2020-6-7
@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
Ameer Hamza 2020-6-3
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
AVM 2020-6-3
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
Ameer Hamza 2020-6-3
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
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
AVM 2020-6-7
@Ameer: Okay the issue is not relevant here. I just asked this question by mistake. It will no longer releted to my case now.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by