کلیدستان

نسخه‌ی کامل: خطا در کد متلب
شما در حال مشاهده نسخه آرشیو هستید. برای مشاهده نسخه کامل کلیک کنید.
سلام دوستان عزیز
سوال والبته مشکلی دارم
چند مدتی هست که درگیر یک برنامه متلب هستم
دائما با error روبرو میشوم، و هر تلاشی کردم نتوانستم در رفعشان موفق شوم!
خواهش میکنم دوستان عزیز بهم کمک کنید
قسمتی از برنامه ام را که نتوانستن اجرا کنم اینجا میگذارم
باتشکر

کد:
a0=1.42;
a1=3/2*a0*px+sqrt(3)/2*a0*py;
a2=3/2*a0*px-sqrt(3)/2*a0*py;
h=6.58479278e-16;
 
kp=k0*exp(-1*(a1*0+a2*0)*i)+k1*exp(-1*(a1*-1+a2*0)*i)+k2*exp(-1*(a1*0+a2*-1)*i)+k3*exp(-1*(a1*1+a2*-1)*i)+k4*exp(-1*(a1*1+a2*0)*i)+k5*exp(-1*(a1*0+a2*1)*i)+k6*exp(-1*(a1*-1+a2*1)*i)+k7*exp(-1*(a1*-1+a2*-1)*i)+k8*exp(-1*(a1*0+a2*-2)*i)+k9*exp(-1*(a1*1+a2*-2)*i)+k10*exp(-1*(a1*-2+a2*1)*i)+k11*exp(-1*(a1*-2+a2*0)*i)+k12*exp(-1*(a1*1+a2*1)*i)+k13*exp(-1*(a1*0+a2*2)*i)
+k14*exp(-1*(a1*-1+a2*2)*i)+k15*exp(-1*(a1*2+a2*-1)*i)+k16*exp(-1*(a1*2+a2*0)*i);
 
avo=6.0221415e23;
mass=-12/avo;
v=3.00e10;
M=[mass  0 0 0 0 0
    0 mass 0 0 0 0
    0 0 mass 0 0 0
    0 0 0 mass 0 0
    0 0 0 0 mass 0
    0 0 0 0 0 mass];
d=gcd(n,m);
if mod((n-m),3*d)==0;
    dr=3*d;
else
    dr=d;
end
a=sqrt(3)*a0;
l=a*sqrt(n^2+m^2+n*m);
dt=1/pi;
rt=dt/2;
t1=(2*m+n)/dr;
t2=-(m+2*n)/dr;
T=sqrt(3)*1/dr;
 
nn=2*(n^2+m^2+n*m)/dr;
b1x=2*pi/(sqrt(3)*a);
b1y=2*pi/a;
b2x=2*pi/(sqrt(3)*a);
b2y=-2*pi/a;
 
q1x=1/nn*(-t2*b1x+t1*b2x);
q1y=1/nn*(-t2*b1y+t1*b2y);
q2x=1/nn*(m*b1x-n*b2x);
q2y=1/nn*(m*b1y-n*b2y);
 
Q2=sqrt(q2x^2+q2y^2);
 
refine=500;
counter=0;
qmin=0;
qmax=pi/T;
qgrid=(qmax-qmin)/(refine-1);
 
for jj=1:refine;
   
    q=qmin+qgrid*(jj-1);
   
    wavevec(jj,1)=q;
   
    w1(jj,1)=q;
    eigv1(jj,1)=q;
   
    w2(jj,1)=q;
    eigv2(jj,1)=q;
   
    w3(jj,1)=q
    eigv3(jj,1)=q;
    
    w4(jj,1)=q
    eigv4(jj,1)=q;
   
    w5(jj,1)=q;
    eigv5(jj,1)=q;
   
    w6(jj,1)=q
    eigv6(jj,1)=q;
end
for ii=0:(nn/2);
    for jj=1:refine
       
q=qmin+qgrid*(jj-1);
 
qx=q/Q2*q2x+ii*q1x;
 
qy=q/Q2*q2y+ii*q1y;
 
kmmm=subs(kp,px,qx);
 
K=subs(kmmm,py,qy);
 
[v,w]=eig(K,M);
 
for mm=1:5
    for mmm=(mm+1):6
       
        if abs(w(mm,mm)) > abs(w(mmm,mmm))
            swap=w(mm,mm);
            w(mm,mm)=w(mmm,mmm);
            w(mmm,mmm)=swap;
           
            for xxx=1:6
                swapv(xxx)=eigv(xxx,mm);
                eigv(xxx,mm)=eigv(xxx,mmm);
                eigv(xxx,mmm)=swapv(xxx);
            end
        end
    end
end
  
for mm=1:6
    sum=0;
    for mmm=1:6
        sum=sum+(abs(eigv(mmm,mm)))^2;
    end
   
    sum=sqrt(sum);
    for mmm=1:6
        eigv(mmm,mm)=eigv(mmm,mm)/sum;
    end
end
 
w1(jj,(ii+2))=real(sqrt(w(1,1)))*h;
for mm=1:6
    eigv1(jj,ii*6+mm+1)=abs(eigv(mm,1));
end
 
w2(jj,(ii+2))=real(sqrt(w(2,2)))*h;
for mm=1:6
    eigv2(jj,ii*6+mm+1)=abs(eigv(mm,2));
end
 
w3(jj,(ii+2))=real(sqrt(w(3,3)))*h;
for mm=1:6
    eigv3(jj,ii*6+mm+1)=abs(eigv(mm,3));
end
 
w4(jj,(ii+2))=real(sqrt(w(4,4)))*h;
for mm=1:6
    eigv4(jj,ii*6+mm+1)=abs(eigv(mm,4));
end
 
w5(jj,(ii+2))=real(sqrt(w(5,5)))*h;
for mm=1:6;
    eigv5(jj,ii*6+mm+1)=abs(eigv(mm,5));
end
 
w6(jj,(ii+2))=real(sqrt(w(6,6)))*h;
for mm=1:6
    eigv6(jj,ii*6+mm+1)=abs(eigv(mm,6));
end
 
counter=counter+1;
percent=counter/(refine*(nn/2+1));
    end
end
 
سلام.
کدهاتون رو کامل اینجا بنویسید. چون من کدهای ارسال قبلیتون رو اجرا کردم و با خطای زیر روبرو شدم :

کد:
??? Undefined function or variable 'px'.

که یعنی متغیر px رو از قبل تعریف نکردید. حدس میزنم بخشی از کد رو یادتون رفته اینجا بنویسید. به هر حال، اگر کدها رو کامل کنید، پس از اجرای کدها، شاید بتونم کمکتون کنم.
با سلام و تشکر
کد کامل را قرار دادم :

کد:
clear all
 clc
 syms theta px py;
 
m=10;
n=10;
 
 
t=[cos(theta) sin(theta) 0
-sin(theta) cos(theta) 0
0 0 1];
 
kk1=[36.5e4 0 0
    0 24.5e4 0
    0 0 9.82e4];
kk2=[8.8e4 0 0
    0 -3.23e4 0
    0 0 -0.4e4];
kk3=[3.000 0 0
    0 -5.25e4 0
    0 0 0.15e4];
kk4=[-1.92e4 0 0
    0 2.29e4 0
    0 0 -0.58e4];
 
 
 
syms kaa theta t;
subs(t,theta,-150);
kaa=t'*kk2*t;
aa=subs(kaa, theta,-150);
aa1=-1*aa;
 
 syms kab theta t;
 subs(t,theta,-120);
kab=t'*kk1*t;
ab=subs(kab,theta,-120);
ab1=-1*ab;
 
 syms kba theta t;
 subs(t,theta,-161);
kba=t'*kk4*t;
ba=subs(kba,theta,-161);
ba1=-1*ba
 
syms kbb theta t;
subs(t,theta,-150)
kbb=t'*kk2*t;
bb=subs(kbb,theta,-150);
bb1=-1*bb;
 
 
k1=[aa ab
ba bb];
 
 syms kaa theta t;
 subs(t,theta,150);
kaa=t'*kk2*t;
aa=subs(kaa,theta,150);
aa2=-1*aa;
 
 syms kab theta t;
 subs(t,theta,120);
kab=t'*kk1*t;
ab=subs(kab,theta,120);
ab2=-1*ab;
 
 syms kba theta t;
 subs(t,theta,161);
kba=t'*kk4*t;
ba=subs(kab,theta,161);
ba2=-1*ba;
 
syms kbb theta t;
subs(t,theta,150);
kbb=t'*kk2*t;
bb=subs(kbb,theta,150);
bb2=-1*bb;
 
k2=[aa ab
    ba bb];
 
 syms kaa theta t;
 subs(t,theta,90);
kaa=t'*kk2*t;
aa=subs(kaa,theta,90);
aa3=-1*aa;
 
 syms kab theta t;
 subs(t,theta,60);
kab=t'*kk3*t;
ab=subs(kab,theta,60);
ab3=-1*ab;
 
 syms kba theta t;
 subs(t,theta,120);
kba=t'*kk3*t;
ba=subs(kba,theta,120);
ba3=-1*ba;
 
syms kbb theta t;
 subs(t,theta,90);
kbb=t'*kk2*t;
subs(kba,theta,90);
bb3=-1*bb;
 
k3=[aa ab
    ba bb];
 
 syms kaa theta t;
 subs(t,theta,30);
kaa=t'*kk2*t;
aa=subs(kaa,theta,30);
aa4=-1*aa;
 
 syms kab theta t;
 subs(t,theta,19);
kab=t'*kk4*t;
ab=subs(kab,theta,19);
ab4=-1*ab;
 
 
 syms kba theta t;
 subs(t,theta,60);
kba=t'*kk1*t;
ba=subs(kba,theta,60);
ba4=-1*ba;
 
 
 syms kbb theta t;
 subs(t,theta,30);
kbb=t'*kk2*t;
bb=subs(kbb,theta,30);
bb4=-1*bb;
 
k4=[aa ab
    ba bb];
 
 
 syms kaa theta t;
 subs(t,theta,-30);
kaa=t'*kk2*t;
aa=subs(kaa,theta,-30);
aa5=-1*aa;
 
syms kab theta t;
 subs(t,theta,-19);
kab=t'*kk4*t;
ab=subs(kab,theta,-19);
ab5=-1*ab;
 
 
 syms kba theta t;
 subs(t,theta,-60);
kba=t'*kk1*t;
ba=subs(kba,theta,-60);
ba5=-1*ba;
 
 syms kbb theta t;
 subs(t,theta,-30);
kbb=t'*kk2*t;
bb=subs(kbb,theta,-30);
bb5=-1*bb;
 
k5=[aa ab
    ba bb];
 
syms kaa theta t;
 subs(t,theta,-90);
kaa=t'*kk2*t;
aa=subs(kaa,theta,-90);
aa6=-1*aa;
 
 syms kab theta t;
 subs(t,theta,-60);
kab=t'*kk3*t;
ab=subs(kba,theta,-60);
ab6=-1*ab;
 
syms kba theta t;
 subs(t,theta,-120);
kba=t'*kk3*t;
ba=subs(kba,theta,-120);
ba6=-1*ba;
 
syms kbb theta t;
 subs(t,theta,-90);
kbb=t'*kk2*t;
bb=subs(kbb,theta,-90);
bb6=-1*bb;
 
k6=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa7=-1*aa;
 
 syms kab theta t;
 subs(t,theta,180);
kab=t'*kk3*t;
ab=subs(kab,theta,180);
ab7=-1*ab;
 
 
ba=[0 0 0
    0 0 0
    0 0 0];
ba7=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb7=-1*bb;
 
 
k7=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa8=-1*aa;
 
 syms kab theta t;
 subs(t,theta,139);
kab=t'*kk4*t;
ab=subs(kab,theta,139);
ab8=-1*ab;
 
 
ba=[0 0 0
    0 0 0
    0 0 0];
ba8=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb8=-1*bb;
 
k8=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa9=-1*aa;
 
 syms kab theta t;
 subs(t,theta,101);
kab=t'*kk4*t;
ab=subs(kab,theta,101);
ab9=-1*ab;
 
 
ba=[0 0 0
    0 0 0
    0 0 0];
ba9=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb9=-1*bb;
 
k9=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa10=-1*aa;
 
 syms kab theta t;
 subs(t,theta,90);
kab=t'*kk4*t;
ab=subs(kab,theta,90);
ab10=-1*ab;
 
 
ba=[0 0 0
    0 0 0
    0 0 0];
ba10=-1*ba
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb10=-1*bb;
 
k10=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa11=-1*aa;
 
 syms kab theta t;
 subs(t,theta,-139);
kab=t'*kk4*t;
ab=subs(kab,theta,-139);
ab11=-1*ab;
 
 
ba=[0 0 0
    0 0 0
    0 0 0];
ba11=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb11=-1*bb;
 
k11=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa12=-1*aa;
 
 
ab=[0 0 0
    0 0 0
    0 0 0];
ab12=-1*ab;
 
 syms kba theta t;
 subs(t,theta,0);
kba=t'*kk3*t;
ba=subs(kba,theta,0);
ba12=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb12=-1*bb;
 
k12=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa13=-1*aa;
 
 
ab=[0 0 0
    0 0 0
    0 0 0];
ab13=-1*ab;
 
 syms kba theta t;
 subs(t,theta,-41);
kba=t'*kk4*t;
ba=subs(kba,theta,-41);
ba13=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb13=-1*bb;
 
 
k13=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa14=-1*aa;
 
 
ab=[0 0 0
    0 0 0
    0 0 0];
ab14=-1*ab;
 
 syms kba theta t;
 subs(t,theta,-79);
kba=t'*kk4*t;
ba=subs(kba,theta,-79);
ba14=-1*ba;
 
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb14=-1*bb;
 
k14=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa15=-1*aa;
 
 
ab=[0 0 0
    0 0 0
    0 0 0];
ab15=-1*ab;
 
 syms kba theta t;
 subs(t,theta,79);
kba=t'*kk4*t;
ba=subs(kba,theta,79);
ba15=-1*ba;
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb15=-1*bb;
 
k15=[aa ab
    ba bb];
 
aa=[0 0 0
    0 0 0
    0 0 0];
aa16=-1*aa;
 
ab=[0 0 0
    0 0 0
    0 0 0];
ab16=-1*ab;
 
 syms kba theta t;
 subs(t,theta,41);
kba=t'*kk4*t;
ba=subs(kba,theta,41);
ba16=-1*ba;
 
bb=[0 0 0
    0 0 0
    0 0 0];
bb16=-1*bb;
 
k16=[aa ab
    ba bb];
 
 
syms kab theta t;
 subs(t,theta,0);
kab=t'*kk1*t;
ab=subs(kab,theta,0);
ab0=-1*ab;
 
 syms kba theta t;
 subs(t,theta,180);
kba=t'*kk1*t;
ba=subs(kba,theta,180);
ba0=-1*ba;
 
aa=aa1+aa2+aa3+aa4+aa5+aa6+aa7+aa8+aa9+aa10+aa11+aa12+aa13+aa14+aa15+aa16+ab0+ab1+ab2+ab3+ab4+ab5+ab6+ab7+ab8+ab9+ab10+ab11+ab12+ab13+ab14+ab15+ab16
bb=ba0+ba1+ba2+ba3+ba4+ba5+ba6+ba7+ba8+ba9+ba10+ba11+ba12+ba13+ba14+ba15+ba16+bb1+bb2+bb3+bb4+bb5+bb6+bb7+bb8+bb9+bb10+bb11+bb12+bb13+bb14+bb15+bb16;
 
k0=[aa ab
    ba bb];
 
a0=1.42;
a1=3/2*a0*px+sqrt(3)/2*a0*py;
a2=3/2*a0*px-sqrt(3)/2*a0*py;
h=6.58479278e-16;
 
kp=k0*exp(-1*(a1*0+a2*0)*i)+k1*exp(-1*(a1*-1+a2*0)*i)+k2*exp(-1*(a1*0+a2*-1)*i)+k3*exp(-1*(a1*1+a2*-1)*i)+k4*exp(-1*(a1*1+a2*0)*i)+k5*exp(-1*(a1*0+a2*1)*i)+k6*exp(-1*(a1*-1+a2*1)*i)+k7*exp(-1*(a1*-1+a2*-1)*i)+k8*exp(-1*(a1*0+a2*-2)*i)+k9*exp(-1*(a1*1+a2*-2)*i)+k10*exp(-1*(a1*-2+a2*1)*i)+k11*exp(-1*(a1*-2+a2*0)*i)+k12*exp(-1*(a1*1+a2*1)*i)+k13*exp(-1*(a1*0+a2*2)*i)
+k14*exp(-1*(a1*-1+a2*2)*i)+k15*exp(-1*(a1*2+a2*-1)*i)+k16*exp(-1*(a1*2+a2*0)*i);
 
avo=6.0221415e23;
mass=-12/avo;
v=3.00e10;
M=[mass  0 0 0 0 0
    0 mass 0 0 0 0
    0 0 mass 0 0 0
    0 0 0 mass 0 0
    0 0 0 0 mass 0
    0 0 0 0 0 mass];
d=gcd(n,m);
if mod((n-m),3*d)==0;
    dr=3*d;
else
    dr=d;
end
a=sqrt(3)*a0;
l=a*sqrt(n^2+m^2+n*m);
dt=1/pi;
rt=dt/2;
t1=(2*m+n)/dr;
t2=-(m+2*n)/dr;
T=sqrt(3)*1/dr;
 
nn=2*(n^2+m^2+n*m)/dr;
b1x=2*pi/(sqrt(3)*a);
b1y=2*pi/a;
b2x=2*pi/(sqrt(3)*a);
b2y=-2*pi/a;
 
q1x=1/nn*(-t2*b1x+t1*b2x);
q1y=1/nn*(-t2*b1y+t1*b2y);
q2x=1/nn*(m*b1x-n*b2x);
q2y=1/nn*(m*b1y-n*b2y);
 
Q2=sqrt(q2x^2+q2y^2);
 
refine=500;
counter=0;
qmin=0;
qmax=pi/T;
qgrid=(qmax-qmin)/(refine-1);
 
for jj=1:refine;
   
    q=qmin+qgrid*(jj-1);
   
    wavevec(jj,1)=q;
   
    w1(jj,1)=q;
    eigv1(jj,1)=q;
   
    w2(jj,1)=q;
    eigv2(jj,1)=q;
   
    w3(jj,1)=q
    eigv3(jj,1)=q;
    
    w4(jj,1)=q
    eigv4(jj,1)=q;
   
    w5(jj,1)=q;
    eigv5(jj,1)=q;
   
    w6(jj,1)=q
    eigv6(jj,1)=q;
end
for ii=0:(nn/2);
    for jj=1:refine
       
q=qmin+qgrid*(jj-1);
 
qx=q/Q2*q2x+ii*q1x;
 
qy=q/Q2*q2y+ii*q1y;
 
kmmm=subs(kp,px,qx);
 
K=subs(kmmm,py,qy);
 
[eigv,w]=eig(K,M);
 
for mm=1:5
    for mmm=(mm+1):6
       
        if abs(w(mm,mm)) > abs(w(mmm,mmm))
            swap=w(mm,mm);
            w(mm,mm)=w(mmm,mmm);
            w(mmm,mmm)=swap;
           
            for xxx=1:6
                swapv(xxx)=eigv(xxx,mm);
                eigv(xxx,mm)=eigv(xxx,mmm);
                eigv(xxx,mmm)=swapv(xxx);
            end
        end
    end
end
  
for mm=1:6
    sum=0;
    for mmm=1:6
        sum=sum+(abs(eigv(mmm,mm)))^2;
    end
   
    sum=sqrt(sum);
    for mmm=1:6
        eigv(mmm,mm)=eigv(mmm,mm)/sum;
    end
end
 
w1(jj,(ii+2))=real(sqrt(w(1,1)))*h;
for mm=1:6
    eigv1(jj,ii*6+mm+1)=abs(eigv(mm,1));
end
 
w2(jj,(ii+2))=real(sqrt(w(2,2)))*h;
for mm=1:6
    eigv2(jj,ii*6+mm+1)=abs(eigv(mm,2));
end
 
w3(jj,(ii+2))=real(sqrt(w(3,3)))*h;
for mm=1:6
    eigv3(jj,ii*6+mm+1)=abs(eigv(mm,3));
end
 
w4(jj,(ii+2))=real(sqrt(w(4,4)))*h;
for mm=1:6
    eigv4(jj,ii*6+mm+1)=abs(eigv(mm,4));
end
 
w5(jj,(ii+2))=real(sqrt(w(5,5)))*h;
for mm=1:6;
    eigv5(jj,ii*6+mm+1)=abs(eigv(mm,5));
end
 
w6(jj,(ii+2))=real(sqrt(w(6,6)))*h;
for mm=1:6
    eigv6(jj,ii*6+mm+1)=abs(eigv(mm,6));
end
 
counter=counter+1;
percent=counter/(refine*(nn/2+1));
    end
end
 
کدها رو کمی تغییر دادم تا دو خطا برطرف شود :

کد:
clear all
  clc
  syms theta px py;
  
m=10;
n=10;
  
  
t=[cos(theta) sin(theta) 0
-sin(theta) cos(theta) 0
0 0 1];
  
kk1=[36.5e4 0 0
     0 24.5e4 0
     0 0 9.82e4];
kk2=[8.8e4 0 0
     0 -3.23e4 0
     0 0 -0.4e4];
kk3=[3.000 0 0
     0 -5.25e4 0
     0 0 0.15e4];
kk4=[-1.92e4 0 0
     0 2.29e4 0
     0 0 -0.58e4];
  
  
  
syms kaa theta t;
subs(t,theta,-150);
kaa=t'*kk2*t;
aa=subs(kaa, theta,-150);
aa1=-1*aa;
  
  syms kab theta t;
  subs(t,theta,-120);
kab=t'*kk1*t;
ab=subs(kab,theta,-120);
ab1=-1*ab;
  
  syms kba theta t;
  subs(t,theta,-161);
kba=t'*kk4*t;
ba=subs(kba,theta,-161);
ba1=-1*ba
  
syms kbb theta t;
subs(t,theta,-150)
kbb=t'*kk2*t;
bb=subs(kbb,theta,-150);
bb1=-1*bb;
  
  
k1=[aa ab
ba bb];
  
  syms kaa theta t;
  subs(t,theta,150);
kaa=t'*kk2*t;
aa=subs(kaa,theta,150);
aa2=-1*aa;
  
  syms kab theta t;
  subs(t,theta,120);
kab=t'*kk1*t;
ab=subs(kab,theta,120);
ab2=-1*ab;
  
  syms kba theta t;
  subs(t,theta,161);
kba=t'*kk4*t;
ba=subs(kab,theta,161);
ba2=-1*ba;
  
syms kbb theta t;
subs(t,theta,150);
kbb=t'*kk2*t;
bb=subs(kbb,theta,150);
bb2=-1*bb;
  
k2=[aa ab
     ba bb];
  
  syms kaa theta t;
  subs(t,theta,90);
kaa=t'*kk2*t;
aa=subs(kaa,theta,90);
aa3=-1*aa;
  
  syms kab theta t;
  subs(t,theta,60);
kab=t'*kk3*t;
ab=subs(kab,theta,60);
ab3=-1*ab;
  
  syms kba theta t;
  subs(t,theta,120);
kba=t'*kk3*t;
ba=subs(kba,theta,120);
ba3=-1*ba;
  
syms kbb theta t;
  subs(t,theta,90);
kbb=t'*kk2*t;
subs(kba,theta,90);
bb3=-1*bb;
  
k3=[aa ab
     ba bb];
  
  syms kaa theta t;
  subs(t,theta,30);
kaa=t'*kk2*t;
aa=subs(kaa,theta,30);
aa4=-1*aa;
  
  syms kab theta t;
  subs(t,theta,19);
kab=t'*kk4*t;
ab=subs(kab,theta,19);
ab4=-1*ab;
  
  
  syms kba theta t;
  subs(t,theta,60);
kba=t'*kk1*t;
ba=subs(kba,theta,60);
ba4=-1*ba;
  
  
  syms kbb theta t;
  subs(t,theta,30);
kbb=t'*kk2*t;
bb=subs(kbb,theta,30);
bb4=-1*bb;
  
k4=[aa ab
     ba bb];
  
  
  syms kaa theta t;
  subs(t,theta,-30);
kaa=t'*kk2*t;
aa=subs(kaa,theta,-30);
aa5=-1*aa;
  
syms kab theta t;
  subs(t,theta,-19);
kab=t'*kk4*t;
ab=subs(kab,theta,-19);
ab5=-1*ab;
  
  
  syms kba theta t;
  subs(t,theta,-60);
kba=t'*kk1*t;
ba=subs(kba,theta,-60);
ba5=-1*ba;
  
  syms kbb theta t;
  subs(t,theta,-30);
kbb=t'*kk2*t;
bb=subs(kbb,theta,-30);
bb5=-1*bb;
  
k5=[aa ab
     ba bb];
  
syms kaa theta t;
  subs(t,theta,-90);
kaa=t'*kk2*t;
aa=subs(kaa,theta,-90);
aa6=-1*aa;
  
  syms kab theta t;
  subs(t,theta,-60);
kab=t'*kk3*t;
ab=subs(kba,theta,-60);
ab6=-1*ab;
  
syms kba theta t;
  subs(t,theta,-120);
kba=t'*kk3*t;
ba=subs(kba,theta,-120);
ba6=-1*ba;
  
syms kbb theta t;
  subs(t,theta,-90);
kbb=t'*kk2*t;
bb=subs(kbb,theta,-90);
bb6=-1*bb;
  
k6=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa7=-1*aa;
  
  syms kab theta t;
  subs(t,theta,180);
kab=t'*kk3*t;
ab=subs(kab,theta,180);
ab7=-1*ab;
  
  
ba=[0 0 0
     0 0 0
     0 0 0];
ba7=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb7=-1*bb;
  
  
k7=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa8=-1*aa;
  
  syms kab theta t;
  subs(t,theta,139);
kab=t'*kk4*t;
ab=subs(kab,theta,139);
ab8=-1*ab;
  
  
ba=[0 0 0
     0 0 0
     0 0 0];
ba8=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb8=-1*bb;
  
k8=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa9=-1*aa;
  
  syms kab theta t;
  subs(t,theta,101);
kab=t'*kk4*t;
ab=subs(kab,theta,101);
ab9=-1*ab;
  
  
ba=[0 0 0
     0 0 0
     0 0 0];
ba9=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb9=-1*bb;
  
k9=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa10=-1*aa;
  
  syms kab theta t;
  subs(t,theta,90);
kab=t'*kk4*t;
ab=subs(kab,theta,90);
ab10=-1*ab;
  
  
ba=[0 0 0
     0 0 0
     0 0 0];
ba10=-1*ba
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb10=-1*bb;
  
k10=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa11=-1*aa;
  
  syms kab theta t;
  subs(t,theta,-139);
kab=t'*kk4*t;
ab=subs(kab,theta,-139);
ab11=-1*ab;
  
  
ba=[0 0 0
     0 0 0
     0 0 0];
ba11=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb11=-1*bb;
  
k11=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa12=-1*aa;
  
  
ab=[0 0 0
     0 0 0
     0 0 0];
ab12=-1*ab;
  
  syms kba theta t;
  subs(t,theta,0);
kba=t'*kk3*t;
ba=subs(kba,theta,0);
ba12=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb12=-1*bb;
  
k12=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa13=-1*aa;
  
  
ab=[0 0 0
     0 0 0
     0 0 0];
ab13=-1*ab;
  
  syms kba theta t;
  subs(t,theta,-41);
kba=t'*kk4*t;
ba=subs(kba,theta,-41);
ba13=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb13=-1*bb;
  
  
k13=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa14=-1*aa;
  
  
ab=[0 0 0
     0 0 0
     0 0 0];
ab14=-1*ab;
  
  syms kba theta t;
  subs(t,theta,-79);
kba=t'*kk4*t;
ba=subs(kba,theta,-79);
ba14=-1*ba;
  
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb14=-1*bb;
  
k14=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa15=-1*aa;
  
  
ab=[0 0 0
     0 0 0
     0 0 0];
ab15=-1*ab;
  
  syms kba theta t;
  subs(t,theta,79);
kba=t'*kk4*t;
ba=subs(kba,theta,79);
ba15=-1*ba;
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb15=-1*bb;
  
k15=[aa ab
     ba bb];
  
aa=[0 0 0
     0 0 0
     0 0 0];
aa16=-1*aa;
  
ab=[0 0 0
     0 0 0
     0 0 0];
ab16=-1*ab;
  
  syms kba theta t;
  subs(t,theta,41);
kba=t'*kk4*t;
ba=subs(kba,theta,41);
ba16=-1*ba;
  
bb=[0 0 0
     0 0 0
     0 0 0];
bb16=-1*bb;
  
k16=[aa ab
     ba bb];
  
  
syms kab theta t;
  subs(t,theta,0);
kab=t'*kk1*t;
ab=subs(kab,theta,0);
ab0=-1*ab;
  
  syms kba theta t;
  subs(t,theta,180);
kba=t'*kk1*t;
ba=subs(kba,theta,180);
ba0=-1*ba;
  
aa=aa1+aa2+aa3+aa4+aa5+aa6+aa7+aa8+aa9+aa10+aa11+aa12+aa13+aa14+aa15+aa16+ab0+ab1+ab2+ab3+ab4+ab5+ab6+ab7+ab8+ab9+ab10+ab11+ab12+ab13+ab14+ab15+ab16
bb=ba0+ba1+ba2+ba3+ba4+ba5+ba6+ba7+ba8+ba9+ba10+ba11+ba12+ba13+ba14+ba15+ba16+bb1+bb2+bb3+bb4+bb5+bb6+bb7+bb8+bb9+bb10+bb11+bb12+bb13+bb14+bb15+bb16;
  
k0=[aa ab
     ba bb];
  
a0=1.42;
a1=3/2*a0*px+sqrt(3)/2*a0*py;
a2=3/2*a0*px-sqrt(3)/2*a0*py;
h=6.58479278e-16;
  
kp=k0*exp(-1*(a1*0+a2*0)*i)+k1*exp(-1*(a1*-1+a2*0)*i)+k2*exp(-1*(a1*0+a2*-1)*i)+k3*exp(-1*(a1*1+a2*-1)*i)+k4*exp(-1*(a1*1+a2*0)*i)+k5*exp(-1*(a1*0+a2*1)*i)+k6*exp(-1*(a1*-1+a2*1)*i)+k7*exp(-1*(a1*-1+a2*-1)*i)+k8*exp(-1*(a1*0+a2*-2)*i)+k9*exp(-1*(a1*1+a2*-2)*i)+k10*exp(-1*(a1*-2+a2*1)*i)+k11*exp(-1*(a1*-2+a2*0)*i)+k12*exp(-1*(a1*1+a2*1)*i)+k13*exp(-1*(a1*0+a2*2)*i)
+k14*exp(-1*(a1*-1+a2*2)*i)+k15*exp(-1*(a1*2+a2*-1)*i)+k16*exp(-1*(a1*2+a2*0)*i);
  
avo=6.0221415e23;
mass=-12/avo;
v=3.00e10;
M=[mass  0 0 0 0 0
     0 mass 0 0 0 0
     0 0 mass 0 0 0
     0 0 0 mass 0 0
     0 0 0 0 mass 0
     0 0 0 0 0 mass];
d=gcd(n,m);
if mod((n-m),3*d)==0;
     dr=3*d;
else
     dr=d;
end
a=sqrt(3)*a0;
l=a*sqrt(n^2+m^2+n*m);
dt=1/pi;
rt=dt/2;
t1=(2*m+n)/dr;
t2=-(m+2*n)/dr;
T=sqrt(3)*1/dr;
  
nn=2*(n^2+m^2+n*m)/dr;
b1x=2*pi/(sqrt(3)*a);
b1y=2*pi/a;
b2x=2*pi/(sqrt(3)*a);
b2y=-2*pi/a;
  
q1x=1/nn*(-t2*b1x+t1*b2x);
q1y=1/nn*(-t2*b1y+t1*b2y);
q2x=1/nn*(m*b1x-n*b2x);
q2y=1/nn*(m*b1y-n*b2y);
  
Q2=sqrt(q2x^2+q2y^2);
  
refine=500;
counter=0;
qmin=0;
qmax=pi/T;
qgrid=(qmax-qmin)/(refine-1);
  
for jj=1:refine;
    
     q=qmin+qgrid*(jj-1);
    
     wavevec(jj,1)=q;
    
     w1(jj,1)=q;
     eigv1(jj,1)=q;
    
     w2(jj,1)=q;
     eigv2(jj,1)=q;
    
     w3(jj,1)=q;
     eigv3(jj,1)=q;
    
     w4(jj,1)=q;
     eigv4(jj,1)=q;
    
     w5(jj,1)=q;
     eigv5(jj,1)=q;
    
     w6(jj,1)=q;
     eigv6(jj,1)=q;
end
for ii=0:(nn/2);
     for jj=1:refine
        
q=qmin+qgrid*(jj-1);
  
qx=q/Q2*q2x+ii*q1x;
  
qy=q/Q2*q2y+ii*q1y;
  
kmmm=subs(kp,px,qx);
  
K=subs(kmmm,py,qy);
  
[eigv,w]=eig(K,M);
  
for mm=1:5
     for mmm=(mm+1):6
        
         if abs(w(mm,mm)) > abs(w(mmm,mmm))
             swap=w(mm,mm);
             w(mm,mm)=w(mmm,mmm);
             w(mmm,mmm)=swap;
            
             for xxx=1:6
                 swapv(xxx)=eigv(xxx,mm);
                 eigv(xxx,mm)=eigv(xxx,mmm);
                 eigv(xxx,mmm)=swapv(xxx);
             end
         end
     end
end
  
for mm=1:6
     sum=0;
     for mmm=1:6
         sum=sum+(abs(eigv(mmm,mm)))^2;
     end
    
     sum=sqrt(sum);
     for mmm=1:6
         eigv(mmm,mm)=eigv(mmm,mm)/sum;
     end
end
  
w1(jj,(ii+2))=real(sqrt(w(1,1)))*h;
for mm=1:6
     eigv1(jj,ii*6+mm+1)=abs(eigv(mm,1));
end
  
w2(jj,(ii+2))=real(sqrt(w(2,2)))*h;
for mm=1:6
     eigv2(jj,ii*6+mm+1)=abs(eigv(mm,2));
end
  
w3(jj,(ii+2))=real(sqrt(w(3,3)))*h;
for mm=1:6
     eigv3(jj,ii*6+mm+1)=abs(eigv(mm,3));
end
  
w4(jj,(ii+2))=real(sqrt(w(4,4)))*h;
for mm=1:6
     eigv4(jj,ii*6+mm+1)=abs(eigv(mm,4));
end
  
w5(jj,(ii+2))=real(sqrt(w(5,5)))*h;
for mm=1:6;
     eigv5(jj,ii*6+mm+1)=abs(eigv(mm,5));
end
  
w6(jj,(ii+2))=real(sqrt(w(6,6)))*h;
for mm=1:6
     eigv6(jj,ii*6+mm+1)=abs(eigv(mm,6));
end
  
counter=counter+1;
percent=counter/(refine*(nn/2+1));
     end
end

اما سپس خطای زیر را می دهد :

کد:
??? Error using ==> sym.eig
Too many input arguments.

Error in ==> Untitled at 560
[eigv,w]=eig(K,M);

وقتی مقادیر  K و M را بررسی می کنیم، نتیجه به صورت زیر است :

کد:
K =

[ -1344806*t*conj(t),             -theta,            -theta,  1344806*t*conj(t),             theta,             theta]
[             -theta, -1012400*t*conj(t),            -theta,              theta, 1012400*t*conj(t),             theta]
[             -theta,             -theta, -361000*t*conj(t),              theta,             theta,  361000*t*conj(t)]
[  1421609*t*conj(t),              theta,             theta, -1364009*t*conj(t),            -theta,            -theta]
[              theta,   868300*t*conj(t),             theta,             -theta, -937000*t*conj(t),            -theta]
[              theta,              theta,  385700*t*conj(t),             -theta,            -theta, -368300*t*conj(t)]


M =

  1.0e-022 *

   -0.1993         0         0         0         0         0
         0   -0.1993         0         0         0         0
         0         0   -0.1993         0         0         0
         0         0         0   -0.1993         0         0
         0         0         0         0   -0.1993         0
         0         0         0         0         0   -0.1993

مشاهده می کنید که M یک ماتریس است اما K یک cell array می باشد (راهنمای دستور cell را بخوانید). بنابراین وقتی دستور eig ، متغیر K را دریافت می کند و می بیند که ماتریس نیست، پیام خطا می دهد. بنابراین باید K را تبدیل به ماتریس بکنید.
ممنون از راهنمایی خوبتون
ولی یه سوال
اون مقادیری که پایینkنوشته شده، شما k را به ماتریس تبدیل کرده اید؟
 
یه فکری به ذهنم زد!
میتونم از ماتریس واحد برای تبدیل k به ماتریس استفاده کنم!?
فقط نمیدانم ماتریس واحد را چطور بنویسم!؟
[1 0 0;0 1 0;0 0 1] اینگونه نوشتم جواب نداد! :
نه من هیچ تغییری در ماتریس K ندادم.
نیازی به این روش ها نیست. تنها باید بدانید که چگونه به مقادیر cell array ارجاع بدهید (چگونه مقادیر را فراخوانی کنید)، سپس مقادیر را از cell array بخوانید و در ماتریسی با نام جدید، ذخیره کنید. آنگاه مقادیر متغیر K را حذف کرده و مقادیر ماتریس جدید را درون آن قرار بدهید.