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;
