function FEM_onedimpbg


%    比較兩種離散方式對精度，計算量的影響。方法1，先代入bloch條件，后離散；方法2，先離散，后代入bloch條件
%u  特向量，bloch函數(shù)的周期部分，第一列表示波數(shù)，第二列表示能帶，第三列表示特征向量
clear
clc

ro=[1;1];
E=[1;1];
d=1;
%參數(shù)
L=2*d;

N=20;

N_band=4;
%N+1個節(jié)點，總共2N+1個
h=d/N;

%實空間單胞數(shù) 2*q

q=10;

x=[0:h:L];
%波數(shù)

mm=20;

b=pi/L/mm;
kk=[-1*pi/L+b/2:b:1*pi/L-b/2];
kk=[kk(1)-b,kk,kk(2*mm)+b];


NN=size(kk);
omiga=[];


ik=E/h;
im=ro*h/6;


u=zeros(NN(2),2*N+1,2*N+1);
%先代入bloch后離散
for l=1:1:NN(2) 
    k=kk(l);
    K=zeros(2*N+1,2*N+1);
    M=zeros(2*N+1,2*N+1);
    
    for j=1:N
        K(j,j)=(2/3*(k*h)^2+2)*ik(1);
        %         K(j,j+1)=(1/6*(k*h)^2+i*k*h-1)*ik(1);
        
        K(j,j+1)=(1/6*(k*h)^2-i*k*h-1)*ik(1);
        
        K(j+1,j)=K(j,j+1)';
        M(j,j)=4*im(1);
        M(j,j+1)=im(1);
        M(j+1,j)=im(1);
    end
    K(1,1)=(1/3*(k*h)^2+1)*ik(1);
    M(1,1)=2*im(1);
    K(N+1,N+1)=(1/3*(k*h)^2+1)*(ik(1)+ik(2));
    K(N,N+1)=K(N+1,N)';
    %     K(N+1,N+2)=(1/6*(k*h)^2+i*k*h-1)*ik(2);
    
    K(N+1,N+2)=(1/6*(k*h)^2-i*k*h-1)*ik(2);
    
    K(N+2,N+1)=K(N+1,N+2)';
    M(N+1,N+1)=2*(im(1)+im(2));
    M(N+1,N+2)=im(2);
    M(N+2,N+1)=im(2);
    for j=N+2:2*N
        K(j,j)=(2/3*(k*h)^2+2)*ik(2);
        %         K(j,j+1)=(1/6*(k*h)^2+i*k*h-1)*ik(2);
        
        K(j,j+1)=(1/6*(k*h)^2-i*k*h-1)*ik(2);
        
        K(j+1,j)=K(j,j+1)';
        M(j,j)=4*im(2);
        M(j,j+1)=im(2);
        M(j+1,j)=im(2);
    end
    K(2*N+1,2*N+1)=(1/3*(k*h)^2+1)*ik(2);
    M(2*N+1,2*N+1)=2*im(2);
    
    KK=K;
    MM=M; 
    
    MM(:,1)=MM(:,1)+MM(:,2*N+1);
    MM(1,:)=MM(1,:)+MM(2*N+1,:);
    KK(:,1)=KK(:,1)+KK(:,2*N+1);
    KK(1,:)=KK(1,:)+KK(2*N+1,:);
    KK(:,2*N+1)=[];
    MM(:,2*N+1)=[];
    KK(2*N+1,:)=[];
    MM(2*N+1,:)=[]; 
    
    
    [eigvetor,eigval]=eig(KK,MM);
    [D1,II]=sort(real(diag(eigval)));
    eigvetor=eigvetor(:,II);
    
    %特征值
    omiga=[omiga,sqrt(D1)];
    
    %特向量
    
    
    %對邊界條件的處理
    eigvetor=[eigvetor;eigvetor(2*N,:)];
    eigvetor=[eigvetor,eigvetor(:,1)];
    
    %歸一化
    for m=1:2*N+1
        eigvetor(:,m)=eigvetor(:,m)/sqrt(eigvetor(:,m)'*eigvetor(:,m));
    end
    
    
    u(l,:,:)=eigvetor;    
end 
omiga=omiga';
y1=[];
for j=1:N_band
    
    y1=[y1,omiga(:,j)];
end 






%先離散
uu=u;
omiga=[];
%先離散,后代入bloch
k=0;
K=zeros(2*N+1,2*N+1);
M=zeros(2*N+1,2*N+1);

for j=1:N
    K(j,j)=(2/3*(k*h)+2)*ik(1);
    K(j,j+1)=(1/6*(k*h)^2+i*k*h-1)*ik(1);
    K(j+1,j)=K(j,j+1)';
    M(j,j)=4*im(1);
    M(j,j+1)=im(1);
    M(j+1,j)=im(1);
end
K(1,1)=(1/3*(k*h)+1)*ik(1);
M(1,1)=2*im(1);
K(N+1,N+1)=(1/3*(k*h)+1)*(ik(1)+ik(2));
K(N,N+1)=K(N+1,N)';
K(N+1,N+2)=(1/6*(k*h)^2+i*k*h-1)*ik(2);
K(N+2,N+1)=K(N+1,N+2)';
M(N+1,N+1)=2*(im(1)+im(2));
M(N+1,N+2)=im(2);
M(N+2,N+1)=im(2);
for j=N+2:2*N
    K(j,j)=(2/3*(k*h)+2)*ik(2);
    K(j,j+1)=(1/6*(k*h)^2+i*k*h-1)*ik(2);
    K(j+1,j)=K(j,j+1)';
    M(j,j)=4*im(2);
    M(j,j+1)=im(2);
    M(j+1,j)=im(2);
end
K(2*N+1,2*N+1)=(1/3*(k*h)+1)*ik(2);
M(2*N+1,2*N+1)=2*im(2);

for l=1:1:NN(2) 
    k=kk(l);
    
    KK=K;
    MM=M; 
    MM(:,1)=MM(:,1)+MM(:,2*N+1)*exp(-i*k*2*d)';
    MM(1,:)=MM(1,:)+MM(2*N+1,:)*exp(-i*k*2*d);
    KK(:,1)=KK(:,1)+KK(:,2*N+1)*exp(-i*k*2*d)';
    KK(1,:)=KK(1,:)+KK(2*N+1,:)*exp(-i*k*2*d);
    
    KK(:,2*N+1)=[];
    MM(:,2*N+1)=[];
    KK(2*N+1,:)=[];
    MM(2*N+1,:)=[]; 
    
    
    [eigvetor,eigval]=eig(KK,MM);
    [D1,II]=sort(real(diag(eigval)));
    eigvetor=eigvetor(:,II);
    omiga=[omiga,sqrt(D1)]; 
    
    eigvetor=[eigvetor;eigvetor(2*N,:)*exp(i*k*L)];
    eigvetor=[eigvetor,eigvetor(:,1)*exp(i*k*L)];
    for m=1:2*N+1
        eigvetor(:,m)=eigvetor(:,m)/sqrt(eigvetor(:,m)'*eigvetor(:,m));
    end
    
    uu(l,:,:)=eigvetor;      
end 
omiga=omiga';
y2=[];
for j=1:6
    
    y2=[y2,omiga(:,j)];
end 

% Uw=wannier(u,kk,b,x,L,N_band);
% Uw=uu;

figure
plot(kk,y1')
hold on
plot(kk,y2','r');