1 function [eigf,eigv,dof]=laplaceeig(node,elem,problem) 2 % 0-boundary eigenvalue problem 3 % problem='0-boundary' 4 [bdNode,Dirichet,inNode]=findboundary(elem); 5 %找出邊界節點、邊界邊的節點編號、內部節點 6 N=size(node,1); %節點個數 7 NT=size(elem,1); %單元個數 8 Nin=N-size(bdNode,1); %內部節點個數 9 ve1=node(elem(:,3),:)-node(elem(:,2),:); 10 %單元第三個節點坐標減去第二個坐標 11 ve2=node(elem(:,1),:)-node(elem(:,3),:); 12 %單元第一個節點坐標減去第三個坐標 13 ve3=node(elem(:,2),:)-node(elem(:,1),:); 14 %單元第二個節點坐標減去第一個坐標 15 area=0.5*(ve3(:,1).*ve2(:,2)+ve3(:,2).*ve2(:,1)) %求每個單元的面積 16 Dlambda(1:NT,:,1)=[-ve1(:,2)./(2*area),ve1(:,1)./(2*area)]; 17 % 求 dlambda1/dx,dlambda1/dy 18 Dlambda(1:NT,:,2)=[-ve2(:,2)./(2*area),ve2(:,1)./(2*area)]; 19 % 求 dlambda2/dx,dlambda2/dy 20 Dlambda(1:NT,:,3)=[-ve3(:,2)./(2*area),ve3(:,1)./(2*area)]; 21 % 求 dlambda3/dx,dlambda3/dy 22 clear ve1,ve2,ve3 23 A=sparse(N,N); %生成N*N的零矩陣 24 B=sparse(N,N); %生成N*N的零矩陣 25 for i=1:3 26 for j=i:3 27 Aij=(Dlambda(:,1,i).*Dlambda(:,1,j)+Dlambda(:,2,i).*Dlambda(:,2,j)).*area; 28 % 求出去 Aij的出去對角的上三角剖分 29 if(j==i) 30 A=A+sparse(elem(:,i),elem(:,j),Aij,N,N);% 求Aij的對角元 31 B=B+sparse(elem(:,i),elem(:,j),area/6,N,N); % lambda^2在邊上的積分 32 else 33 A=A+sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],[Aij;Aij],N,N); 34 B=B+sparse([elem(:,i);elem(:,j)],[elem(:,j);elem(:,i)],[area;area]/12,N,N); 35 %lambda_i*lambda_j在單元的積分 36 end 37 end 38 39 end 40 clear Aij 41 % 0-boundary 注意節點個數不夠需要進行一致加密 42 A(bdNode,:)=[];A(:,bdNode)=[]; 43 B(bdNode,:)=[];B(:,bdNode)=[]; %去掉邊界節點 44 [eigf,eigv]=eigs(A,B,1,'sm'); %按模最小求第一個特征值及對應的特征向量 45 eigf=eigf/(eigf'*A*eigf)^0.5; %按1模把特征值向量規范化 u/(u'*A*u) 46 eigf=accumarray([inNode,ones(Nin,1)],eigf,[N,1]);% 表 u 組裝成與節點個數大小 47 dof=Nin;