Eigenvalues of Elliptical Operators over Arbitrary Regions

Richard Gass

Department of Physics University of Cincinnati

Abstract
Mathematica 10 introduced a number of powerful functions for region based computation. NDSolve can now solve PDE’s over arbitrary regions but currently built-in methods cannot solve elliptical eigenvalue problems over arbitrary regions. I will show how to use the new region and mesh functions to solve elliptical eigenvalue problems over arbitrary regions. These methods are very fast and efficient and allow the solution of quite complex problems. Several physically important examples will be discussed.

◀     |     ▶

Solving PDE’s over regions

In Version 10 it is easy to solve a PDE over a region

ℛ=RegionDifference[Disk[{0,0},1],Disk[{0,0},1/4]];
f=NDSolveValue[{(Subscript[∇, {x,y}]^2)u[x,y]==0,DirichletCondition[u[x,y]==2,(x^2)+(y^2)==1],DirichletCondition[u[x,y]==3,(x^2)+(y^2)==1/16]},u,{x,y}∈ℛ];
GraphicsGrid[{{RegionPlot[ℛ,PlotLabel->"The region ℛ",AxesLabel->{"x","y"}],Plot3D[f[x,y],{x,y}∈ℛ,PlotLabel->"Solution of (∇^2)U=0",AxesLabel->{"x","y"}]}},ImageSize->6 72]

Built-in methods can solve an eigenvalue problem over a region. However the new geometric mesh functions make this relatively easy using finite element methods

◀     |     ▶

Eigenvalue problems over regions

We want to solve (∇^2)ψ+(k^2)ψ=0 over a region ℛ.

• First generate the weak form of the equation: (Underscript[∫∫, ℛ])(∇ (Subscript[T, s]) · ∇ ψ) ds=(k^2)(Underscript[∫∫, ℛ])(Subscript[T, s])  ψ where (Subscript[T, s]) is a test function.
• Discretize ℛ with first order triangular elements

• On each element ψ =(∑_(i=1)^3)(Subscript[α, i])(x,y) (Subscript[ψ, i])  where (Subscript[ψ, i])=(Subscript[a, i])+(Subscript[b , i])x+(Subscript[c, i])y with (Subscript[a, i])=(Subscript[x, j])(Subscript[y, k])-(Subscript[x, k])(Subscript[y, j]), (Subscript[b, i])=(Subscript[y, j])-(Subscript[y, k]) , (Subscript[c, i])=(Subscript[x, k])-(Subscript[x, j]) and (Subscript[α, i])(x,y)=(1/(2A))((Subscript[a, i])+(Subscript[b, i]) x+(Subscript[c, i])y) with A the area of the triangular element.

• Pick (Subscript[T, s])=α(x,y) . This is the usual Galerkin method.

◀     |     ▶

• Result is the for each element we have (Subscript[S, △])ψ=(k^2)(Subscript[T, △])ψ with (Subscript[S, △])=(Underscript[∫∫, △])∇(Subscript[α, i])·∇(Subscript[α, j])dxdy and (Subscript[T, △])= (Underscript[∫∫, △])(Subscript[α, i])·(Subscript[α, j])dx dy where the subscript △ indicates the triangular element. (Subscript[S, △]) and (Subscript[T, △]) are three by three matrices

• Assemble the local matrices (Subscript[S, △]) and (Subscript[T, △]) into global n ×n matrices S and T where n is the total number of nodes

• Solve the generalized eigenvalue problem S ψ =(k^2) T ψ.

◀     |     ▶

Example - a hexagon

Generate and triangulate a  boundary mesh region

ngon[p_,q_]:=Table[{Cos[2Pi k q/p],Sin[2Pi k q/p]},{k,p}]
boundaryMeshRegion=BoundaryMeshRegion[ngon[6,1],Line[{1,2,3,4,5,6,1}]]

ℳ=TriangulateMesh[boundaryMeshRegion,MeshQualityGoal->"Maximal"]

We will also need the boundary in order to impose the boundary conditions

ℬ=BoundaryMesh[ℳ];

◀     |     ▶

Legended[Show[HighlightMesh[ℳ,Labeled[0,"Index"]],Graphics[{Black,PointSize[0.01],Map[Point,MeshCoordinates[ℬ]]}],Graphics[{Red,Map[Point,DeleteDuplicates[Delete[MeshCoordinates[ℳ],Flatten[Map[Position[MeshCoordinates[ℳ],#]&,MeshCoordinates[ℬ]],1]]]]}],ImageSize->7 72,PlotLabel->"Vertices"],Column[{Row[{Graphics[{{Black,PointSize[0.1],Point[{0,0}]},Style[Text["Boundary Vertex",{3/4,0}],18]}]}],Row[{Graphics[{{Red,PointSize[0.1],Point[{0,0}]},Style[Text["Free Vertex",{3/4,0}],18]}]}]}]]

◀     |     ▶

We need the coordinates of the vertices of each triangle which we can get using MeshCoordinates and MeshCells.

coords=Partition[Map[Part[MeshCoordinates[ℳ],#]&,Flatten[Map[Apply[List,#]&,MeshCells[ℳ,2]]]],3];

We can now evaluate the element functions on each triangle and compute (Subscript[S, △]) and (Subscript[T, △]) . The matrix (Subscript[T, △]) depends only on the area of the triangle.

b=Map[{#[[2,2]]-#[[2,3]],#[[2,3]]-#[[2,1]],#[[2,1]]-#[[2,2]]}&,Map[Transpose,coords]];
c=Map[{#[[1,3]]-#[[1,2]],#[[1,1]]-#[[1,3]],#[[1,2]]-#[[1,1]]}&,Map[Transpose,coords]];
a=Map[Cross[#1[[1]],#1[[2]]]&,Map[Transpose,coords]];
A=(1/2)Map[Det,coords/.{x_,y_}->{1,x,y}];
element[b_,c_]:={{((b[[1]])^2)+((c[[1]])^2),b[[1]]b[[2]]+c[[1]]c[[2]],b[[1]]b[[3]]+c[[1]]c[[3]]},{b[[2]]b[[1]]+c[[2]]c[[1]],((b[[2]])^2)+((c[[2]])^2),b[[2]]b[[3]]+c[[2]]c[[3]]},{b[[3]]b[[1]]+c[[3]]c[[1]],b[[3]]b[[2]]+c[[3]]c[[2]],((b[[3]])^2)+((c[[3]])^2)}};
Sel=1/(4 A) MapThread[element[#1,#2]&,{b,c}];
Tel=Map[Times[{{2,1,1},{1,2,1},{1,1,2}},#]&,A/12];

◀     |     ▶

Assemble the global matrix

Best understood with a very small example.

𝒥=MeshRegion[{{0., 0.},{0.5,0},{1.0,0.0} ,{0.75,0.5},{0.25,0.5},{0.56,1}},{
    Polygon[{{1,2,5},{2,3,4},{2,4,5},{5,4,6}}]}];
HighlightMesh[𝒥,{Style[Labeled[0,"Index"],Red],Labeled[2,"Index"]}]

We can extract the mesh cells and the indexes of the mesh cells. The indices are the vertices of the triangles in counterclockwise order.

MeshCells[𝒥,2]

indexes=MeshCells[𝒥,2]/.Polygon[x_]->x

We can generate the list of global indices by taking the outer product of the indices with it’s self.

globalIndexList=Map[Outer[List,#,#]&,indexes]

Dimensions[globalIndexList]

For each element (triangle) we have a 3 x 3  matrix. To assemble the global matrix we map (Subscript[S, ij]^e)-> (Subscript[S, n(i,e),n(j,e)]) where e denotes the element number and n(e,i) is the global node number indexed by the element number e and the local node number. Thus n(2,3)=5 since that is the number of the second node (going counterclockwise) of the third element. This is done with the function generateMatrix.

generateMatrix[index_List,elementMatrix_List]:=SparseArray[Flatten[MapThread[Rule,{index,elementMatrix},2]],dimension,0];
dimension=Length[MeshCoordinates[𝒥]];
S=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalIndexList,Sel}]];
T=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalIndexList,Tel}]];

◀     |     ▶

Bundling the code up

The function drop (taken from taken from stack exchange) is a generalization of Drop . It is used to drop the nodes that are fixed by the boundary conditions.

drop[m_,parts__List]/;Length@{parts}≤ArrayDepth[m]:=m[[##]]&@@MapThread[Complement,{Range@Dimensions[m,Length@{parts}],{parts}}]

t1=TimeUsed[];coords=Partition[Map[Part[MeshCoordinates[ℳ],#]&,Flatten[Map[Apply[List,#]&,MeshCells[ℳ,2]]]],3];
b=Map[{#[[2,2]]-#[[2,3]],#[[2,3]]-#[[2,1]],#[[2,1]]-#[[2,2]]}&,Map[Transpose,coords]];
c=Map[{#[[1,3]]-#[[1,2]],#[[1,1]]-#[[1,3]],#[[1,2]]-#[[1,1]]}&,Map[Transpose,coords]];
a=Map[Cross[#1[[1]],#1[[2]]]&,Map[Transpose,coords]];
element[b_,c_]:={{((b[[1]])^2)+((c[[1]])^2),b[[1]]b[[2]]+c[[1]]c[[2]],b[[1]]b[[3]]+c[[1]]c[[3]]},{b[[2]]b[[1]]+c[[2]]c[[1]],((b[[2]])^2)+((c[[2]])^2),b[[2]]b[[3]]+c[[2]]c[[3]]},{b[[3]]b[[1]]+c[[3]]c[[1]],b[[3]]b[[2]]+c[[3]]c[[2]],((b[[3]])^2)+((c[[3]])^2)}};
A=(1/2)Map[Det,coords/.{x_,y_}->{1,x,y}];
Sel=1/(4 A) MapThread[element[#1,#2]&,{b,c}];
Tel=Map[Times[{{2,1,1},{1,2,1},{1,1,2}},#]&,A/12];
indexes=MeshCells[ℳ,2]/.Polygon[x_]->x;
globalIndexList=Map[Outer[List,#,#]&,indexes];
dimension=Length[MeshCoordinates[ℳ]];
generateMatrix[index_List,elementMatrix_List]:=SparseArray[Flatten[MapThread[Rule,{index,elementMatrix},2]],dimension,0];
S=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalIndexList,Sel}]];
T=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalIndexList,Tel}]];
boundaryNodes=Flatten[Map[Position[MeshCoordinates[ℳ],#]&,MeshCoordinates[ℬ]]];
boundaryPoints=MeshCoordinates[ℬ];
𝒮=drop[S,boundaryNodes,boundaryNodes];
𝒯=drop[T,boundaryNodes,boundaryNodes];
bc=0.0;
bcData=Transpose[{DeleteDuplicates[boundaryPoints],Table[bc,{i,1,Length[boundaryNodes]}]}];
meshCoords=MeshCoordinates[ℳ];
freePoints=DeleteDuplicates[Delete[meshCoords,Flatten[Map[Position[meshCoords,#]&,boundaryPoints],1]]];
solutions=Sort[Transpose[Eigensystem[{𝒮,𝒯}]]];
freeData=Transpose[{freePoints,Take[solutions[[1,2]],{1,Length[freePoints]}]}];
TimeUsed[]-t1

data=Partition[Flatten[Join[bcData,freeData]],3];
ListDensityPlot[data/.{x_,y_,z_}->{x,y,(z^2)},PlotLabel->"Ground State"]

◀     |     ▶

A more complicated example: From a picture

Picture

First we import the picture. We will specify a potential on the inner hexagon.

picture=Import["/Users/gass/Desktop/HexWithWell.JPG"];

◀     |     ▶

Then we select points on the outer and inner boundaries.

outerBoundaryData={{583.2517985611511`,353.66906474820166`},{593.6115107913669`,330.8776978417268`},{597.7553956834532`,312.23021582733827`},{601.8992805755396`,301.8705035971225`},{610.1870503597122`,287.3669064748203`},{622.6187050359712`,274.93525179856135`},{630.9064748201439`,254.21582733812966`},{647.4820143884892`,252.14388489208648`},{664.0575539568346`,252.14388489208648`},{697.2086330935251`,254.21582733812966`},{734.5035971223022`,252.14388489208648`},{763.5107913669065`,245.92805755395705`},{802.8776978417266`,239.7122302158275`},{838.1007194244604`,239.7122302158275`},{856.7482014388489`,241.78417266187068`},{871.2517985611511`,254.21582733812966`},{887.8273381294963`,272.8633093525182`},{902.3309352517985`,289.4388489208635`},{918.9064748201439`,308.086330935252`},{933.410071942446`,330.8776978417268`},{952.0575539568346`,349.5251798561153`},{966.5611510791367`,370.2446043165469`},{956.2014388489209`,384.74820143884904`},{945.8417266187051`,399.2517985611513`},{931.338129496403`,428.2589928057556`},{914.7625899280575`,461.41007194244617`},{896.115107913669`,488.3453237410074`},{875.3956834532374`,519.4244604316549`},{862.9640287769785`,540.1438848920865`},{848.4604316546762`,556.7194244604318`},{833.956834532374`,562.9352517985612`},{817.3812949640288`,567.0791366906476`},{792.5179856115108`,567.0791366906476`},{771.7985611510792`,567.0791366906476`},{738.6474820143885`,575.3669064748203`},{715.8561151079136`,575.3669064748203`},{688.9208633093525`,577.4388489208635`},{659.9136690647482`,575.3669064748203`},{647.4820143884892`,573.2949640287771`},{637.1223021582733`,565.0071942446045`},{624.6906474820144`,556.7194244604318`},{608.115107913669`,540.1438848920865`},{591.5395683453237`,515.2805755395685`},{574.9640287769785`,494.56115107913683`},{558.3884892086331`,469.6978417266189`},{543.884892086331`,448.9784172661872`},{550.1007194244604`,440.6906474820146`},{560.4604316546762`,424.11510791366925`},{570.8201438848921`,403.39568345323755`},{574.9640287769785`,376.4604316546764`}};

innerBoundaryData={{711.7122302158274`,382.67625899280586`},{728.2877697841726`,378.5323741007195`},{742.7913669064748`,378.5323741007195`},{765.5827338129496`,376.4604316546763`},{792.5179856115108`,374.38848920863313`},{817.3812949640288`,372.31654676258995`},{829.8129496402878`,382.67625899280586`},{838.1007194244604`,399.2517985611511`},{852.6043165467626`,422.04316546762595`},{865.0359712230215`,446.9064748201439`},{846.3884892086331`,457.26618705035975`},{825.6690647482014`,469.6978417266188`},{796.6618705035971`,486.27338129496405`},{778.0143884892086`,496.6330935251799`},{761.4388489208633`,496.6330935251799`},{751.0791366906475`,488.34532374100723`},{736.5755395683453`,475.91366906474826`},{722.0719424460432`,465.5539568345324`},{701.3525179856115`,453.12230215827344`},{697.2086330935251`,430.33093525179856`},{693.0647482014389`,407.5395683453238`},{695.136690647482`,386.8201438848921`}};

We need to make the lists of points cyclical so that they define boundary regions.

cyclicOuterBoundaryData=Join[outerBoundaryData,{First[outerBoundaryData]}];
cyclicInnerBoundaryData=Join[innerBoundaryData,{First[innerBoundaryData]}];
ListPlot[{cyclicOuterBoundaryData,cyclicInnerBoundaryData},AspectRatio->1,PlotRange->{{500,1000},{200,600}},Axes->{True,True}]

We now construct the two boundaries

ℬ1=BoundaryMeshRegion[cyclicOuterBoundaryData,Line[Join[Range[1,Length[outerBoundaryData]],{1}]]];
ℬ2=BoundaryMeshRegion[cyclicInnerBoundaryData,Line[Join[Range[1,Length[innerBoundaryData]],{1}]]];

Scale the data

The next step is to scale the boundaries. This is done by calculating a scaling function that scales the data so that the x coordinates lie between xmin and xmax.

Clear[scalingFunction,f]
scalingFunction[d_List,{xStart_,xEnd_}]:=Module[{xmin,xmax},
xmin=Min[d/.{x_,y_}->x];
xmax=Max[d/.{x_,y_}->x];
First[(c1+c2 x)/.Solve[{c1+c2 xmin==xStart,c1+c2 xmax==xEnd},{c1,c2}]]]
f[x_]=scalingFunction[cyclicOuterBoundaryData,{-1,1}];

We now scale both boundaries with the same scaling function.

ℬ1=BoundaryMeshRegion[Partition[Map[f,Flatten[cyclicOuterBoundaryData]],2],Line[Join[Range[1,Length[outerBoundaryData]],{1}]]];
ℬ2=BoundaryMeshRegion[Partition[Map[f,Flatten[cyclicInnerBoundaryData]],2],Line[Join[Range[1,Length[innerBoundaryData]],{1}]]];

and triangulate the outer boundary.

ℋ3=TriangulateMesh[ℬ1]

◀     |     ▶

Pick out the mesh in the inner region

We need to be able to select mesh cells that lie inside the inner boundary.

The first step is to generate a region function

rm=RegionMember[ℬ2]

Pick out the cells whose centroid lies inside the inner boundary. We get the centroids by using the function PropertyValue. For example:

PropertyValue[{ℋ3,2},MeshCellCentroid][[64]]

The cells in the inner region are given by

cellsInRegion=Flatten[Map[Take[MeshCells[ℋ3,2],#]&,Position[rm[PropertyValue[{ℋ3,2},MeshCellCentroid]],True]]];

This is a list of all the polygons that lie inside ℬ2. The numbers are the global indices of the vertices.

Take[cellsInRegion,{1,10}]

The cell coordinates for the these cells are given by

cells=Map[Polygon,Partition[Map[Part[MeshCoordinates[ℋ3],#]&,Flatten[Map[Apply[List,#]&,cellsInRegion]]],3]];

We will need these to compute the area of the selected cells. This is a list of polygons but now the polygons are specified by the coordinates of their vertices instead of by the global indices of the vertices.

Take[cells,10]

We can visualize the polygons with

Graphics[cells]

and combine them with the outer boundary region.

Show[ℋ3,Graphics[{Red,Opacity[0.5],cells}]]

This does not do a very good job of capturing the inner boundary because our mesh is too coarse.

◀     |     ▶

Show[ℋ3,Graphics[{Red,Opacity[0.5],cells}],Graphics[Line[Partition[Map[f,Flatten[cyclicInnerBoundaryData]],2]]]]

We need to refine the mesh but we want to do this only around the inner boundary.

meshRefinementFunction=Function[{vertices,area},Block[{x,y},{x,y}=Mean[vertices];If[RegionDistance[ℬ2,{x,y}]<.05&&RegionDistance[RegionDifference[ℬ1,ℬ2],{x,y}]<.05,area>0.0005,area>0.005]]];

ℳ=TriangulateMesh[ℬ1, MeshRefinementFunction->meshRefinementFunction]

The refined mesh does a much better job

cellsInRegion=Flatten[Map[Take[MeshCells[ℳ,2],#]&,Position[rm[PropertyValue[{ℳ,2},MeshCellCentroid]],True]]];
cells=Map[Polygon,Partition[Map[Part[MeshCoordinates[ℳ],#]&,Flatten[Map[Apply[List,#]&,cellsInRegion]]],3]];
Show[ℳ,Graphics[{Red,Opacity[0.25],cells}],Graphics[{Red,Map[Point,Partition[Map[f,Flatten[cyclicOuterBoundaryData]],2]]}],Graphics[{PointSize[.01],Green,Map[Point,Partition[Map[f,Flatten[cyclicInnerBoundaryData]],2]]}],Graphics[Line[Partition[Map[f,Flatten[cyclicInnerBoundaryData]],2]]],ImageSize->5 72]

◀     |     ▶

We will need the global indices of the selected cells and these are given by

dimension=Length[MeshCoordinates[ℳ]];

selectedMeshIndices=cellsInRegion/.Polygon[x_]->x;
globalSelectedIndexList=Map[Outer[List,#,#]&,selectedMeshIndices];

The coordinates of the selected cells are

selectedCells=cells/.Polygon[x_]->x;

areaSelectedCells=(1/2)Map[Det,selectedCells/.{x_,y_}->{1,x,y}];

Since the potential is piecewise constant the matrix elements are given by (Underscript[(Subscript[V, 0])∫∫, Δ])(Subscript[α, i])(Subscript[α, j])dxdy. Note that up to a constant this just (Subscript[T, el]). Thus we can define a matrix (Subscript[V, el]) as

Vel=V0 Map[Times[{{2,1,1},{1,2,1},{1,1,2}},#]&,areaSelectedCells/12];

and assemble it into a global matrix.

V=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalSelectedIndexList,Vel}]];

Set the potential to be attractive

V0=-20;

◀     |     ▶

Code Body

indexes=MeshCells[ℳ,2]/.Polygon[x_]->x;
globalIndexList=Map[Outer[List,#,#]&,indexes];
generateMatrix[index_List,elementMatrix_List]:=SparseArray[Flatten[MapThread[Rule,{index,elementMatrix},2]],dimension,0];
ℬ=BoundaryMesh[ℳ];
coords=Partition[Map[Part[MeshCoordinates[ℳ],#]&,Flatten[Map[Apply[List,#]&,MeshCells[ℳ,2]]]],3];
b=Map[{#[[2,2]]-#[[2,3]],#[[2,3]]-#[[2,1]],#[[2,1]]-#[[2,2]]}&,Map[Transpose,coords]];
c=Map[{#[[1,3]]-#[[1,2]],#[[1,1]]-#[[1,3]],#[[1,2]]-#[[1,1]]}&,Map[Transpose,coords]];
a=Map[Cross[#1[[1]],#1[[2]]]&,Map[Transpose,coords]];
A=(1/2)Map[Det,coords/.{x_,y_}->{1,x,y}];
element[b_,c_]:={{((b[[1]])^2)+((c[[1]])^2),b[[1]]b[[2]]+c[[1]]c[[2]],b[[1]]b[[3]]+c[[1]]c[[3]]},{b[[2]]b[[1]]+c[[2]]c[[1]],((b[[2]])^2)+((c[[2]])^2),b[[2]]b[[3]]+c[[2]]c[[3]]},{b[[3]]b[[1]]+c[[3]]c[[1]],b[[3]]b[[2]]+c[[3]]c[[2]],((b[[3]])^2)+((c[[3]])^2)}};
Sel=1/(4 A) MapThread[element[#1,#2]&,{b,c}];
Tel=Map[Times[{{2,1,1},{1,2,1},{1,1,2}},#]&,A/12];


S=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalIndexList,Sel}]];
T=Apply[Plus,MapThread[generateMatrix[#1,#2]&,{globalIndexList,Tel}]];
boundaryNodes=Flatten[Map[Position[MeshCoordinates[ℳ],#]&,MeshCoordinates[ℬ]]];
boundaryPoints=MeshCoordinates[ℬ];
𝒮=drop[S+V,boundaryNodes,boundaryNodes];
𝒯=drop[T,boundaryNodes,boundaryNodes];
bc=0.0;
bcData=Transpose[{DeleteDuplicates[boundaryPoints],Table[bc,{i,1,Length[boundaryNodes]}]}];
meshCoords=MeshCoordinates[ℳ];
freePoints=DeleteDuplicates[Delete[meshCoords,Flatten[Map[Position[meshCoords,#]&,boundaryPoints],1]]];
solutions=Sort[Transpose[Eigensystem[{𝒮,𝒯}]]];
freeData=Map[Transpose[{freePoints,Take[solutions[[#,2]],{1,Length[freePoints]}]}]&,Range[Length[freePoints]]];
data=Map[Partition[Flatten[Join[bcData,freeData[[#]]]],3]&,Range[Length[freePoints]]];

◀     |     ▶

labels={"Ground State","First Excited State","Second Excited State","Third Excited State","Fourth Excited State"};
Manipulate[Show[ListDensityPlot[data[[i]]/.{x_,y_,z_}->{x,y,(z^2)},PlotLabel->labels[[i]]],Graphics[Line[Partition[Map[f,Flatten[cyclicInnerBoundaryData]],2]]]],{i,1,5,1},SaveDefinitions->True]

◀     |     ▶