版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜

返回列表

菲斐

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 21
帖子: 2
在线: 1.5小时
虫号: 3247358
注册: 2014-06-01
专业: 凝聚态物性 II ：电子结构

[求助] 求matlab中 diff Difference order N must be a positive integer scalar什么意思已有1人参与

我是matlab新手，不知道程序中为什么会出现??? Error using ==> diff
Difference order N must be a positive integer scalar.

Error in ==> chern at 90
phpx2=diff(H2,'x');
我只是把第一个程序中的H变了一下，就有错误，可是如果把这2段分到2个m文件里，就木有问题，咋回事呀。大神帮忙呀，程序如下

syms x y
a=1.0;
t=1;
t1=0.1*t;
u0=0.6*sqrt(3)*t;
u2=0.6*sqrt(3)*t;
f=2*cos(x*a/2)*cos(sqrt(3)*a*y/2)-2*i*cos(sqrt(3)*a*y/2)*sin(a*x/2)+cos(x*a)+i*sin(x*a);
g=-i*(2*i*sin(sqrt(3)*a*y)-4*i*sin(sqrt(3)*a*y/2+3*a*x/2));
H=[u0/2+t1*g-u2/2,2*t*f;2*t*conj(f),-u0/2-t1*g+u2/2];
phpx=diff(H,'x');
phpy=diff(H,'y');
h=6.63*10^(-34);%此为普朗克常数

hall=0;

for kx = -2*pi/(3*a) : 4*pi/(30*a) : 2*pi/(3*a)

for ky = -4*pi*sqrt(3)/(9*a) : 8*pi*sqrt(3)/(90*a) : 4*pi*sqrt(3)/(9*a)

      if -2*pi/(3*a) <= kx <= 0 && (-sqrt(3)*kx/3-4*pi*sqrt(3)/(9*a)) <= ky <=(sqrt(3)*kx/3+4*pi*sqrt(3)/(9*a))

         % f=exp(-i*((kx*a/2)+sqrt(3)*ky*a/2))+ exp(-i*((kx*a/2)-sqrt(3)*ky*a/2))+exp(-i*(-kx*a)) ;
      %g1=exp(i*(ky*sqrt(3)*a))-exp(i*(-kx*3*a/2+ky*sqrt(3)*a/2))+exp(i*(-kx*3*a/2-ky*sqrt(3)*a/2));
      %g2=-exp(i*(-ky*sqrt(3)*a))+exp(i*(kx*3*a /2-ky*sqrt(3)*a/2))-exp(i*(kx*3*a/2+ky*sqrt(3)*a/2));
      %g=-i*(g1+g2);


         f=2*cos(kx*a/2)*cos(sqrt(3)*a*ky/2)-2*i*cos(sqrt(3)*a*ky/2)*sin(a*kx/2)+cos(kx*a)+i*sin(kx*a);
         g=-i*(2*i*sin(sqrt(3)*a*ky)-4*i*sin(sqrt(3)*a*ky/2+3*a*kx/2));
         H=[u0/2+t1*g-u2/2,2*t*f;2*t*conj(f),-u0/2-t1*g+u2/2];
         [P,E]=eig(H);
         mk=[P(1,1);P(2,1)]; %用eig求得的本征矢和本征值不是按从小到大的顺序排列的
         nk=[P(1,2);P(2,2)];
         Em=E(1,1);
         En=E(2,2);
         vx=subs(phpx,{x,y},{kx,ky});
         vy=subs(phpy,{x,y},{kx,ky});
         s=i*exp(2)*9*a^2/(4*sqrt(3)*pi*h);%此为系数
         if Em<En
         ha=s*((mk'*vx*nk)*(nk'*vy*mk)-(mk'*vy*nk)*(nk'*vx*mk)) /(Em-En)^2;

         else
         ha=s*((nk'*vx*mk)*(mk'*vy*nk)-(nk'*vy*mk)*(mk'*vx*nk)) /(En-Em)^2;

         end
         hall=hall+ha;
      end
      if 0 <= kx <= 2*pi/(3*a) && (sqrt(3)*kx/3-4*pi*sqrt(3)/(9*a)) <= ky <=(-sqrt(3)*kx/3+4*pi*sqrt(3)/(9*a))

         %f=exp(-i*((kx*a/2)+sqrt(3)*ky*a/2))+ exp(-i*((kx*a/2)-sqrt(3)*ky*a/2))+exp(-i*(-kx*a)) ;
      %g1=exp(i*(ky*sqrt(3)*a))-exp(i*(-kx*3*a/2+ky*sqrt(3)*a/2))+exp(i*(-kx*3*a/2-ky*sqrt(3)*a/2));
      %g2=-exp(i*(-ky*sqrt(3)*a))+exp(i*(kx*3*a /2-ky*sqrt(3)*a/2))-exp(i*(kx*3*a/2+ky*sqrt(3)*a/2));
      %g=-i*(g1+g2);



         f=2*cos(kx*a/2)*cos(sqrt(3)*a*ky/2)-2*i*cos(sqrt(3)*a*ky/2)*sin(a*kx/2)+cos(kx*a)+i*sin(kx*a);
         g=-i*(2*i*sin(sqrt(3)*a*ky)-4*i*sin(sqrt(3)*a*ky/2+3*a*kx/2));
         H=[u0/2+t1*g-u2/2,2*t*f;2*t*conj(f),-u0/2-t1*g+u2/2];
         [P,E]=eig(H);
         mk=[P(1,1);P(2,1)];
         nk=[P(1,2);P(2,2)];
         Em=E(1,1);
         En=E(2,2);
         vx=subs(phpx,{x,y},{kx,ky});
         vy=subs(phpy,{x,y},{kx,ky});
         s=i*exp(2)*9*a^2/(4*sqrt(3)*pi*h);
         if Em<En
            ha=s*((mk'*vx*nk)*(nk'*vy*mk)-(mk'*vy*nk)*(nk'*vx*mk)) /(Em-En)^2;

         else
            ha=s*((nk'*vx*mk)*(mk'*vy*nk)-(nk'*vy*mk)*(mk'*vx*nk)) /(En-Em)^2;

         end
         hall=hall+ha;


      end
end
end
c1=hall*h/(2*pi*exp(2));%第一段

H2=[u0/2-t1*g+u2/2,2*t*f;2*t*conj(f),-u0/2+t1*g-u2/2];
phpx2=diff(H2,'x');
phpy2=diff(H2,'y');
hall2=0;
for kx = -2*pi/(3*a) : 4*pi/(30*a) : 2*pi/(3*a)

for ky = -4*pi*sqrt(3)/(9*a) : 8*pi*sqrt(3)/(90*a) : 4*pi*sqrt(3)/(9*a)

      if -2*pi/(3*a) <= kx <= 0 && (-sqrt(3)*kx/3-4*pi*sqrt(3)/(9*a)) <= ky <=(sqrt(3)*kx/3+4*pi*sqrt(3)/(9*a))

         % f=exp(-i*((kx*a/2)+sqrt(3)*ky*a/2))+ exp(-i*((kx*a/2)-sqrt(3)*ky*a/2))+exp(-i*(-kx*a)) ;
      %g1=exp(i*(ky*sqrt(3)*a))-exp(i*(-kx*3*a/2+ky*sqrt(3)*a/2))+exp(i*(-kx*3*a/2-ky*sqrt(3)*a/2));
      %g2=-exp(i*(-ky*sqrt(3)*a))+exp(i*(kx*3*a /2-ky*sqrt(3)*a/2))-exp(i*(kx*3*a/2+ky*sqrt(3)*a/2));
      %g=-i*(g1+g2);


         f=2*cos(kx*a/2)*cos(sqrt(3)*a*ky/2)-2*i*cos(sqrt(3)*a*ky/2)*sin(a*kx/2)+cos(kx*a)+i*sin(kx*a);
         g=-i*(2*i*sin(sqrt(3)*a*ky)-4*i*sin(sqrt(3)*a*ky/2+3*a*kx/2));
         H2=[u0/2-t1*g+u2/2,2*t*f;2*t*conj(f),-u0/2+t1*g-u2/2];
         [P,E]=eig(H2);
         mk=[P(1,1);P(2,1)]; %用eig求得的本征矢和本征值不是按从小到大的顺序排列的
         nk=[P(1,2);P(2,2)];
         Em=E(1,1);
         En=E(2,2);
         vx2=subs(phpx2,{x,y},{kx,ky});
         vy2=subs(phpy2,{x,y},{kx,ky});
         s=i*exp(2)*9*a^2/(4*sqrt(3)*pi*h);%此为系数
         if Em<En
         ha=s*((mk'*vxb*nk)*(nk'*vyb*mk)-(mk'*vyb*nk)*(nk'*vxb*mk)) /(Em-En)^2;

         else
         ha=s*((nk'*vxb*mk)*(mk'*vyb*nk)-(nk'*vyb*mk)*(mk'*vxb*nk)) /(En-Em)^2;

         end
         hall2=hallb+ha;
      end
      if 0 <= kx <= 2*pi/(3*a) && (sqrt(3)*kx/3-4*pi*sqrt(3)/(9*a)) <= ky <=(-sqrt(3)*kx/3+4*pi*sqrt(3)/(9*a))

         %f=exp(-i*((kx*a/2)+sqrt(3)*ky*a/2))+ exp(-i*((kx*a/2)-sqrt(3)*ky*a/2))+exp(-i*(-kx*a)) ;
      %g1=exp(i*(ky*sqrt(3)*a))-exp(i*(-kx*3*a/2+ky*sqrt(3)*a/2))+exp(i*(-kx*3*a/2-ky*sqrt(3)*a/2));
      %g2=-exp(i*(-ky*sqrt(3)*a))+exp(i*(kx*3*a /2-ky*sqrt(3)*a/2))-exp(i*(kx*3*a/2+ky*sqrt(3)*a/2));
      %g=-i*(g1+g2);



         f=2*cos(kx*a/2)*cos(sqrt(3)*a*ky/2)-2*i*cos(sqrt(3)*a*ky/2)*sin(a*kx/2)+cos(kx*a)+i*sin(kx*a);
         g=-i*(2*i*sin(sqrt(3)*a*ky)-4*i*sin(sqrt(3)*a*ky/2+3*a*kx/2));
         H2=[u0/2-t1*g+u2/2,2*t*f;2*t*conj(f),-u0/2+t1*g-u2/2];
         [P,E]=eig(H2);
         mk=[P(1,1);P(2,1)];
         nk=[P(1,2);P(2,2)];
         Em=E(1,1);
         En=E(2,2);
         vx2=subs(phpx2,{x,y},{kx,ky});
         vy2=subs(phpy2,{x,y},{kx,ky});
         s=i*exp(2)*9*a^2/(4*sqrt(3)*pi*h);
         if Em<En
            ha=s*((mk'*vxb*nk)*(nk'*vyb*mk)-(mk'*vyb*nk)*(nk'*vxb*mk)) /(Em-En)^2;

         else
            ha=s*((nk'*vxb*mk)*(mk'*vyb*nk)-(nk'*vyb*mk)*(mk'*vxb*nk)) /(En-Em)^2;

         end
         hall2=hall2+ha;


      end
end
end
c2=hallb*h/(2*pi*exp(2));
c=c1+c2;%第二段

回复此楼

» 猜你喜欢

有没有人能给点建议已经有5人回复
假如你的研究生提出不合理要求已经有12人回复
实验室接单子已经有7人回复
全日制（定向）博士已经有5人回复
萌生出自己或许不适合搞科研的想法，现在跑or等等看？已经有4人回复
Materials Today Chemistry审稿周期已经有4人回复
参与限项已经有3人回复
对氯苯硼酸纯化已经有3人回复
所感已经有4人回复
要不要辞职读博？已经有7人回复

1楼 2014-10-30 17:36:52

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zhchh008

金虫 (正式写手)

应助: 13 (小学生)
金币: 1613.2
散金: 861
红花: 13
帖子: 833
在线: 202.6小时
虫号: 99456
注册: 2005-11-11
性别: GG
专业: 资源化工

【答案】应助回帖

感谢参与，应助指数 +1

第一： u0=0.6*sqrt(3)*t; u2=0.6*sqrt(3)*t; 二者是一样的。第二，这是典型的VB或C语言习惯。完全没有matrix的影子。错误信息的意思是，你在求导（微分）时，阶数必须为正整数，比如不能求某函数的0,5阶导数。

赞一下

回复此楼

2楼2014-10-30 19:53:58

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

菲斐

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 21
帖子: 2
在线: 1.5小时
虫号: 3247358
注册: 2014-06-01
专业: 凝聚态物性 II ：电子结构

引用回帖:

2楼: Originally posted by zhchh008 at 2014-10-30 19:53:58
第一： u0=0.6*sqrt(3)*t; u2=0.6*sqrt(3)*t; 二者是一样的。第二，这是典型的VB或C语言习惯。完全没有matrix的影子。错误信息的意思是，你在求导（微分）时，阶数必须为正整数，比如不能求某函数的0,5阶导数。

u0与u2的比值确实就是一样的呀，不一样的只是H和H2,我才学matlab，见谅，我这是想求chern number，您听说过吗？

赞一下

回复此楼

3楼2014-11-04 16:03:16

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

相关版块跳转我要订阅楼主菲斐的主题更新

返回列表

24小时热门版块排行榜

菲斐

[求助] 求matlab中 diff Difference order N must be a positive integer scalar什么意思 已有1人参与

» 猜你喜欢

zhchh008

【答案】应助回帖

菲斐

[求助] 求matlab中 diff Difference order N must be a positive integer scalar什么意思已有1人参与