MATLAB中文论坛

 找回密码
 注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

查看: 291|回复: 1|关注: 0

[未答复] 双重积分出现两个问题,非有限结果和已达计算最大数目

[复制链接]

新手

10 麦片

财富积分


050


3

主题

4

帖子

0

最佳答案
发表于 2017-10-16 20:39:58 | 显示全部楼层 |阅读模式
我要积的双重积分出现了两个问题,一是计算速度很慢,二是警告:“非有限结果。积分未成功。可能具有奇异性”,以及“已达到函数计算的最大数目(10000)。这个结果将使全局误差测试失败”
我的编程用的方法也很入门,请各位大侠指点。

在这里先表感谢!

我要计算的方程是:
                              

QQ截图20171016203230.png
                              
这里的被积函数已经完全编在下面了,里面的G,E,f都是k和theta_k的函数,

function y=ImintegrandW(k,thk)   %这个方程为被积函数

global q thq VF  s1 s2 lambda kB T u eta omega   %要用到的全局参数,eta是方程里面的\delta是一个无穷小量,我赋值为0.05*omega。omega为频率

KQ=sqrt(k.^2+q.^2+2.*k.*q.*cos(thk-thq));          %矢量的模|k+q|

Snkq=(k.*sin(thk)+q.*sin(thq))./KQ;                      %sin(theta_k+q)
Cskq=(k.*cos(thk)+q.*cos(thq))./KQ;                    %cos(theta_k+q)
Cs3thkq=4.*(Cskq).^3-3.*Cskq;                            %cos(3*theta_k+q)

Ekq=s2.*sqrt(VF.^2.*KQ.^2+lambda.^2.*KQ.^6.*Cs3thkq.^2);
Ek=s1.*sqrt(VF.^2.*k.^2+lambda.^2.*k.^6.*(cos(3.*thk)).^2);

Rk=sqrt(VF.^2.*k.^2+lambda.^2.*k.^6.*(cos(3.*thk)).^2);
Phik=atan(sqrt((Rk+s1.*lambda.*k.^3.*(cos(3.*thk)))./(Rk-s1.*lambda.*k.^3.*(cos(3.*thk)))));

Rkq=sqrt(VF.^2.*KQ.^2+lambda.^2.*KQ.^6.*(Cs3thkq).^2);
Phikq=atan(sqrt((Rkq+s2.*lambda.*KQ.^3.*(Cs3thkq))./(Rkq-s2.*lambda.*KQ.^3.*(Cs3thkq))));

Coskq=cos(thk).*Cskq+sin(thk).*Snkq;

%%% overlap factor
Overlap=(sin(Phikq)).^2.*(sin(Phik)).^2+(cos(Phikq)).^2.*(cos(Phik)).^2....
+2.*s1.*s2.*(sin(Phikq).*cos(Phikq)).*(sin(Phik).*cos(Phik)).*Coskq;

%%% delta function
DeltaF=eta./((Ekq-Ek-omega).^2+eta.^2);

%%% Fermi Dirac distribution
Fkq=1./(exp((Ekq-u)./T./kB)+1);
Fk=1./(exp((Ek-u)./T./kB)+1);
Fermi=Fkq-Fk;

y=Overlap.*DeltaF.*Fermi.*k;


主程序如下:
clear; clc;

global q thq VF  s1 s2 lambda kB T u omega eta   %参数全局

kB=8.617e-5;                          % eVK-1;
VF=0.255;                             % 费米速度 eV*nm
lambda=0;                              %参数 eV*nm^3  
u=0.25;                                    %费米能级
T=0;                                         %温度T=0

wR=0:0.01*u:0.01*u;  L_w=length(wR);             %频率段

qR=linspace(0,1,100);              %不同的q的值作为画图的横坐标


thq=0;              %矢量q的角度,设置为0.

s1=-1;  s2=1;  
for i_w=L_w:L_w;     omega=wR(i_w);      eta=0.01*omega;          %本来是做omega的循环的,这里只取一个点做测试

    ImPwth=zeros(1,L_q);

    for i_q=1:L_q; q=qR(i_q);   

              ImPwth(i_q)=integral2(@(k,thk) ImintegrandW(k,thk),0,20,0,2*pi);   %对k和thk(矢量k的角度)积分

    end  %iq

    plot(qR,ImPwth); hold on


end; %%% L_w     



新手

10 麦片

财富积分


050


3

主题

4

帖子

0

最佳答案
 楼主| 发表于 2017-10-16 20:46:02 | 显示全部楼层
计算结果如下: QQ截图20171016204410.png
QQ截图20171016204426.png

您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /1 下一条

快速回复 返回顶部 返回列表