This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 1

Consider the 2D eigenvalue problem -Δ (Subscript[u, λ])(x,y)=λ (Subscript[u, λ])(x,y) on a square S. The boundary conditions on two edges and the two adjacent half edges are homogeneous Dirichlet conditions, meaning  (Subscript[u, λ])(x,y)(Subscript[|, (∂_D)S])=0 (so-called mixed or hybrid boundary conditions [✶], [✶]). On the remaining two half edges that meet at a vertex, the boundary condition is a homogeneous Robin condition (1-α)∂(Subscript[u, λ])(x,y)/∂n+α (Subscript[u, λ])(x,y)(Subscript[|, (∂_D)S])=0. Make an animation showing how some of the higher states change as α varies from 1 to 0 (meaning from homogeneous Dirichlet conditions of the half edge to homogeneous Neumann conditions). Use a finite element discretization to solve the 2D eigenvalue problem.

(∫_S)∇u(x,y).∇v (x,y)dx dy+(α/(1-α))(∫_((∂_D)S)) u(x,y) v(x,y) ds=(Underscript[λ∫, S]) u(x,y) v(x,y) dx dy,

where we assume that the homogeneous Dirichlet conditions are strictly fulfilled.

A convenient discretization of the square is the subdivision into smaller squares. And a convenient shape function basis set are biquadratic polynomials with knot points at the vertices, the edge midpoints, and the center of a square.

Evaluate[Table[p[k], {k, 9}]] = 
           {{0, 0}, {2, 0}, {2, 2}, {0, 2}, {1, 0}, {2, 1}, {1, 2}, {0, 1}, {1, 1}}/2;

Here are the explicit forms of the nine shape functions.

Table[shapeFunction[m, {ξ_, η_}] =
 Factor[Sum[c[i, j] ξ^i η^j, {i, 0, 2}, {j, 0, 2}] /.
   Solve[Table[(Sum[c[i, j] ξ^i η^j, {i, 0, 2}, {j, 0, 2}] /. 
      Thread[{ξ, η} -> p[k]]) == KroneckerDelta[m, k], {k, 9}],  
         Flatten[Table[c[i, j], {i, 0, 2}, {j, 0, 2}]]][[1]]], {m, 9}]

Here are the nine shape functions graphically hown.

Show[GraphicsArray[#]]& /@ Partition[
Table[Plot3D[Evaluate[shapeFunction[m, {ξ, η}]], {ξ, 0, 1}, {η, 0, 1}, 
             PlotRange -> All, ColorFunction -> GrayLevel, Axes -> True,
             DisplayFunction -> Identity], {m, 9}], 3];

Next, we must calculate the mass matrix, the stiffness matrix, and the boundary line integral matrix for our nine shape functions. Later, we will use subsquares of edge length δ, so the Jacobian is (δ^2).

𝒥 = δ^2;

Here is the mass matrix.

ℳℳ = 𝒥 Table[Integrate[shapeFunction[m, {ξ, η}] shapeFunction[n, {ξ, η}],
                         {ξ, 0, 1}, {η, 0, 1}], {m, 9}, {n, 9}];

The calculation of the stiffness matrix is straightforward too.

(* unit square coordinates as a function of real coordinates *)
ξ[{x_, y_}] := (x - x1)/δ (* with δx = δ = (x2 - x1) *);
η[{x_, y_}] := (y - y1)/δ (* with δy = δ = (y4 - y1) *);

(* shape parameter values *)
{𝒜, 𝒞, ℬ} = 
{D[ξ[{x, y}], x]^2 + D[ξ[{x, y}], y]^2,
 D[η[{x, y}], x]^2 + D[η[{x, y}], y]^2,
 D[ξ[{x, y}], x] D[η[{x, y}], x] + D[ξ[{x, y}], y] D[η[{x, y}], y]} 𝒥

𝒮ℳ = 
With[{sf = shapeFunction[#, {ξ, η}]&},
 Table[𝒜 Integrate[D[sf[m], ξ] D[sf[n], ξ], {ξ, 0, 1}, {η, 0, 1}] +
       𝒞 Integrate[D[sf[m], η] D[sf[n], η], {ξ, 0, 1}, {η, 0, 1}],
        {m, 9}, {n, 9}]];

And here is the matrix for the boundary line integrals.

ℬℳ = Table[0, {m, 9}, {n, 9}];

Function[{ξη, args, nodes}, (ℬℳ[[##]] = 
 δ Integrate[shapeFunction[#1, args] shapeFunction[#2, args],
                         {ξη, 0, 1}])& @@@ 
  (* points that are connected *)
  ({{#1, #1}, {#2, #2}, {#3, #3}, {#1, #2}, {#2, #1}, 
    {#1, #3}, {#3, #1}, {#2, #3}, {#3, #2}}& @@ nodes)] @@@
    {{ξ, {ξ, 1}, {3, 4, 7}}, {ξ, {ξ, 0}, {1, 5, 2}},
     {η, {1, η}, {2, 3, 6}}, {η, {0, η}, {1, 4, 8}}};

Next, we subdivide our original square [-1,1]×[-1,1] into subsquares of edge length δ. The function makeSubSquares generates the subsquares and the function pointList maps the knot points into the subsquares.

makeSubSquares[δ_] := 
 Table[Polygon[{{x, y}, {x + δ, y}, {x + δ, y + δ}, {x, y + δ}}],
        {x, -1, 1 - δ, δ}, {y, -1, 1 - δ, δ}] // Flatten[#, 1]&;

pointList[Polygon[{{x1_, y1_}, {x2_, y2_}, {x3_, y3_}, {x4_, y4_}}]] :=
             ({x1, y1} + {#1 (x2 - x1), #2 (y4 - y1)})& @@@ Array[p, 9]

Now, we classify all discretization points into interior, boundary, and Robin-boundary points. The function makeVertexAssignments subdivides the initial square with Robin boundary conditions according to the test robinPointQ (each boundary point is subjected to this test) into (𝒹^2) squares and generates list of vertices for the interior points, and the Robin-boundary points. For later use, it makes various variable names globally visible. This allows to reuse them easily in various functions.

makeVertexAssignments[𝒹_, robinPointQ_] := 
Module[{subSquarePoints, allVertices, boundaryVertices},
(* global variables *)
Clear[subSquares, allNonDirichletVertices, dirichletVertices, 
      robinVertices, λ, δ, nonDirichletVertexQ, robinVertexQ];
δ = 2/𝒹;
(* form subsquares and mapped points *)
subSquares = makeSubSquares[δ];
subSquarePoints = pointList /@ subSquares;
allVertices = Union[Level[subSquarePoints, {-2}]];
(* all boundary vertices *)
boundaryVertices = Cases[allVertices, {-1, _} | {+1, _} | {_, -1} | {_, +1}];
robinVertices = Select[boundaryVertices, robinPointQ];
(* make functions for Robin-boundary points *)
Do[robinVertexQ[robinVertices[[j]]] = True,
   {j, Length[robinVertices]}];
dirichletVertices = Complement[boundaryVertices, robinVertices];
allNonDirichletVertices = Complement[allVertices, dirichletVertices];
λ = Length[allNonDirichletVertices];
(* make functions for inner and Robin-boundary points *)
Do[nonDirichletVertexQ[allNonDirichletVertices[[j]]] = True;
   coordinateVertex[j] = allNonDirichletVertices[[j]];
   vertexNumber[allNonDirichletVertices[[j]]] = j, 
   {j, λ}];
Do[nonDirichletVertexQ[dirichletVertices[[j]]] = False;
   coordinateVertex[λ + j] = dirichletVertices[[j]];
   vertexNumber[dirichletVertices[[j]]] = λ + j, 
   {j, Length[dirichletVertices]}]]

Here is graphic showing the Dirichlet boundary points in green, the Robin boundary points in red, and all non-Dirichlet points in gray. The subsquares are bounded by blue lines.

(* function to visualize point assignments *)
discretizationGraphics[] := 
Show[Graphics[{
 {RGBColor[0, 0, 1], subSquares /. Polygon[l_] :> Line[Append[l, First[l]]]},
 {PointSize[0.02], RGBColor[0, 1, 0], Point /@ dirichletVertices},
 {PointSize[0.012], RGBColor[1, 0, 0], Point /@ robinVertices},
 {PointSize[0.006], GrayLevel[0.4], Point /@ allNonDirichletVertices}}], 
      PlotRange -> All, AspectRatio -> Automatic, Frame -> True];

robinPointQ[{x_, y_}] := (0 < x <= 1 && y == 1) || 
                         (0 < y <= 1 && x == 1)
makeVertexAssignments[12, robinPointQ]

discretizationGraphics[]

For the solution of the eigenvalue problem through finite elements, we must assemble the global mass and stiffness matrices. The function makeSparseMatrix generates the global mass matrix, stiffness matrix, and the boundary integral matrix arising from the Robin boundary conditions. To assemble the matrix, we use the global values of the functions subSquares, vertexNumber, nonDirichletVertexQ, and robinVertexQ. Because the resulting matrices are sparse and large, we generate them as sparse arrays.

(* add elements *)
makeMatrix[values_] := ((List @@ #[[1, 1]]) -> (Plus @@ (Last /@ #)))& /@ 
                        Split[Sort[values], #1[[1]] == #2[[1]]&]

makeSparseMatrix[criteria_, ℳ_] :=
Module[{L = {}},
 (* loop over all subsquares *)
 Do[xyPoints = pointList[subSquares[[k]]];
   (* loop over all nine points in a subsquare *)
   Do[(* for interior and Robin-boundary points, add contributions *)
      If[criteria[xyPoints[[a]]] && criteria[xyPoints[[b]]],
         {i, j} = vertexNumber /@ xyPoints[[{a, b}]];
         L = {L, C[i, j] -> ℳ[[a, b]]}], {a, 9}, {b, 9}],
     {k, Length[subSquares]}];
SparseArray[makeMatrix @ Flatten[L], λ {1, 1}, 0]]

We unite the generation of the three matrices in the function sparseMatrixTriple.

sparseMatrixTriple[] := makeSparseMatrix @@@
 {{nonDirichletVertexQ, 𝒮ℳ}, {nonDirichletVertexQ, ℳℳ}, {robinVertexQ, ℬℳ}}

We obtain matrices of size 552×552 for the above discretization example.

Dimensions /@ sparseMatrixTriple[]

Next, we define some functions for the visualizations of the eigenfunctions. The function eigenFunctionListPlot3D generates a 3D plot of the eigenfunction.

(* z-coordinate for the 3D plot is eigenvector value *)
zCoordinate[{x_, y_}, evec_] := 
            If[nonDirichletVertexQ[{x, y}], evec[[vertexNumber[{x, y}]]], 0]

(* array of z-coordinates *)
valueTable[evec_] := 
Module[{subsubSquarePoints = Table[{x, y}, {x, -1, 1, δ}, {y, -1, 1, δ}]},
       Map[Append[#, zCoordinate[#, evec]]&, subsubSquarePoints, {-2}]]

eigenFunctionListPlot3D[evec_] :=
Module[{subsubSquarePolygons = 
        Table[Polygon[{{x, y}, {x + δ, y}, {x + δ, y + δ}, {x, y + δ}}],
        {x, -1, 1 - δ, δ}, {y, -1, 1 - δ, δ}] // Flatten[#, 1]&, polys},
 (* form four polygons for each subsquare *)
 polys = Map[{#[[1]], #[[2]], zCoordinate[#, evec]}&, 
             subsubSquarePolygons, {-2}];
  (* make 3D graphics *)
 Graphics3D[{EdgeForm[], polys}, 
            PlotRange -> All, BoxRatios -> {1, 1, 0.4}]]

For higher, more rapidly changing states, a contour plot is typically better suited for visualizations. This time, we show (by default) the absolute value of the eigenfunctions through gray levels and the nodal lines as thin red lines.

(* values for cl homogeneously spaced contours *)
homContours[data_, cl_] := 
With[{λ = Length[Flatten[data]]},
    #[[Round[λ/cl]]]& /@ Partition[Sort[Flatten[data]], Round[λ/cl]]]

eigenFunctionListContourPlot[evec_, f_:Abs] :=
Module[{subsubSquarePoints = Table[{x, y}, {y, -1, 1, δ}, {x, -1, 1, δ}],
        data1, data2, lcp1, lcp2},
 (* prepare data *)
 data1 =  Map[zCoordinate[#, evec]&, subsubSquarePoints, {-2}];
 data2 =  Map[f, data1, {2}];
 (* make graphics *)
 Block[{$DisplayFunction = Identity}, 
 data1a = Take[data1, {2, -2}, {2, -2}];
 (* nodal lines *)
 lcp1 = ListContourPlot[data1a, PlotRange -> All, ContourLines -> True, 
                 ContourShading -> False, Contours -> {0},
                 ContourStyle -> {Hue[0], Thickness[0.004]},
                 MeshRange -> {{-1 + δ, 1 - δ}, {-1 + δ, 1 - δ}}];
 (* absolute value (or f of) the eigenfunctions *)
 lcp2 = ListContourPlot[data2, PlotRange -> All, Frame -> True, 
                 ContourLines -> False, MeshRange -> {{-1, 1}, {-1, 1}},
                 Contours -> homContours[data2, Round[Sqrt[Length[evec]]]],
                 ContourStyle -> {Thickness[0.002]}];
 Show[{lcp2, lcp1, (* border *)
       Graphics[{Thickness[0.01], RGBColor[0, 0, 1], 
                  Line[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}]}]},
      Frame -> False]]]

Before implementing the animation, we visualize some example eigenfunctions for various parts of the square’s boundary having Robin boundary conditions. The function exampleBoundaryGraphics generates the mass, stiffness, and boundary matrices, and solves the eigenvalue for the list of α-values αs for the lowest eigenvalues and for eigenvalues of approximate size shift. The lowest eigenfunctions we display as 3D plots, and the higher ones as contour plots.

exampleBoundaryGraphics[robinCondition_, 𝒹_, αs_, shift_] :=
Module[{nEvs = 3, evOpts, A, B, Ar, evals, evecs},
 evOpts = Sequence["Arnoldi", BasisSize -> 50, MaxIterations -> 10^4,
                              Tolerance -> 10^-8];
 (* generate matrices *)
 makeVertexAssignments[𝒹, robinCondition]; 
 {A, B, Ar} = sparseMatrixTriple[];
 (* show problem discretized structure *)
 Show[Block[{$DisplayFunction = Identity}, discretizationGraphics[]],
      Frame -> False];
 (* show higher eigenfunctions for various α-values *)
 Function[α, Print["α = ", N @ α]; ABPair = N[{A + α/(1 - α) Ar, B}];
  (* calculate and show lowest eigenfunctions for various α-values *)
  {evals, evecs} = Reverse /@ 
  Eigensystem[ABPair, 3, Method -> {evOpts, Shift -> 0}];
  Show[GraphicsArray[eigenFunctionListPlot3D /@ evecs]];
  (* calculate and show higher eigenfunctions for various α-values *)
  {evals, evecs} = Reverse /@ 
  Eigensystem[ABPair, nEvs, Method -> {evOpts, Shift -> shift}];
  Show[GraphicsArray[#]]& /@ 
      Partition[eigenFunctionListContourPlot[Abs[#]]& /@ evecs, 3]] /@ αs]

In the first example, the middle third of one edge has Robin boundary conditions.

robinPointQ[{x_, y_}] := y == 1 && -1/3 < x < 1/3;
exampleBoundaryGraphics[robinPointQ, 50, {0, 0.5, 0.99}, 200]

In the second example, one half of each edges has Robin boundary conditions. This time, we show higher states.

robinPointQ[{x_, y_}] := x == -1 && y < 0 || x == +1 && y > 0 ||
                         y == -1 && x > 0 || y == +1 && x < 0;
exampleBoundaryGraphics[robinPointQ, 50, {0, 0.5, 0.99}, 600]

In the third example, two diagonally opposite half edge has Robin boundary conditions.

robinPointQ[{x_, y_}] := x == -1 && y < 0 || x == +1 && y > 0;
exampleBoundaryGraphics[robinPointQ, 50, {0, 0.5, 0.99}, 600]

In the fourth example, the middle thirds of all four edges have Robin boundary conditions.

robinPointQ[{x_, y_}] := 
     x == -1 && -1/3 < y < 1/3 || x == +1 && -1/3 < y < 1/3 ||
     y == -1 && -1/3 < x < 1/3 || y == +1 && -1/3 < x < 1/3;
exampleBoundaryGraphics[robinPointQ, 50, {0, 0.5, 0.99}, 600]

A typical feature of all of the last examples are the existence of relatively large areas of uniform patterns [✶].

And in the last example, we assign the Robin boundary condition randomly to all vertices.

robinPointQ[{x_, y_}] := randomRobinValue[{x, y}];
(* remember which points are Robin-boundary points *)
randomRobinValue[{x_, y_}] := 
      randomRobinValue[{x, y}] = Random[Integer] === 0
exampleBoundaryGraphics[robinPointQ, 50, {0, 0.5, 0.99}, 600]

Now, we return to our animation problem. To get an idea where and how fast the transition from Dirichlet to Neumann boundary condition is reflected in the eigensystems, we use a relatively crude discretization and calculate and visualize the lowest eigenvalues. We parametrize α through α=1-(10^(Subscript[α, E])) to easily span a large range of α-values.

robinPointQ[{x_, y_}] := (0 < x <= 1 && y == 1) || 
                         (0 < y <= 1 && x == 1)
makeVertexAssignments[7, robinPointQ];
{A, B, Ar} = sparseMatrixTriple[];
Dimensions[A]

Show[Graphics[{PointSize[0.003],
Table[α = 1 - 10^αExp;
     Point[{αExp, #}]& /@ 
     (Take[Eigenvalues[N @ Normal[{A + α/(1 - α) Ar, B}]], -12]/Pi^2),  
      {αExp, 0, -2.5, -1/200}]}], Frame -> True, PlotRange -> {-1, 5}]

We see that most of the transition happens for -2<(Subscript[α, E])<0. For Neumann boundary conditions (α=0), all states are nondegenerate, for Dirichlet boundary conditions (α->1), states become degenerate. We also see level crossing (this is possible because it is a 2D problem).

For the actual animation we will use (60^2) subsquares. This results in matrices of size 14280×14280.

robinPointQ[{x_, y_}] := (0 < x <= 1 && y == 1) || 
                         (0 < y <= 1 && x == 1)
makeVertexAssignments[60, robinPointQ];
{A, B, Ar} = sparseMatrixTriple[];
Dimensions[A]

Here are the points visualized.

discretizationGraphics[]

We will use states with energies (λ-values) around 1010. To make sure that the calculated eigensystems are reliable, we compare the calculate spectrum for the nearly-Dirichlet case of α=0.999999999. The eigenvalues deviate about 0.1% from the exact values.

α = 1 - 10.^-9;
evals = Eigenvalues[N @ {A + α/(1 - α) Ar, B}, 12, Method -> {"Arnoldi",   
                                   BasisSize -> 240, MaxIterations -> 10^4,
                                   Tolerance -> 12, Shift -> 1010}]/(Pi^2/4) // Reverse

Select[Sort[Flatten[Table[k^2 + l^2, {k, 50}, {l, 50}]]], (Min[%] - 3 < # < Max[%] + 3)&]

Now, we can implement the animation for varying α. For each α, we calculate a few eigensystems. Because the eigenvalues of the states change with α, we use a dynamically averaged value for the Shift option of Eigensystem. To select the “corresponding” eigenvectors, we cannot rely on the eigenvalues because they become degenerate. So, we form the overlap between the last eigenvector and the newly calculated eigenvectors for the next value of α. The state with the highest overlap is selected. We will show the α-dependence for three states.

We start with the Neumann case.

(* number of eigensystems to calculate for each value of α *) o = 25;
ABPair = {A, B};
eigensystem = Eigensystem[ABPair // N, o, Method -> {"Arnoldi",   
                                   BasisSize -> 240, MaxIterations -> 10^4,
                                   Tolerance -> 10^-8, Shift -> 1010}];
                                   
currentEvals = eigensystem[[1]][[{14, 13, 12}]];
currentEvecs = eigensystem[[2]][[{14, 13, 12}]];

Next, we loop over the α-values and store the resulting contour plots in the list animationFrames.

animationFrames[framesToShow_] := 
With[{allFrames = 60},
Module[{l1 = Table[αExp, {αExp, 0, 3, 1/allFrames}], l2,
        αExpValuesToShow, overlaps, maxOverlapPositions},
(* which α-values to display *)
αExpValuesToShow = If[framesToShow === All, l1, 
    l2 = Table[αExp, {αExp, 0, 3, 3/framesToShow}];
    Function[s, Sort[{N[Abs[s - #]], #}& /@ l1][[1, 2]]] /@ l2]; 
Do[α = 1 - 10^-αExp;
   (* calculate eigenvalues and eigenvectors *)
   ABPair = {A + α/(1 - α) Ar, B};
   currentShift = Plus @@ currentEvals/Length[currentEvals];
   eigensystem = Eigensystem[ABPair // N, o, Method -> {"Arnoldi",   
                                 BasisSize -> 240, MaxIterations -> 10^4,
                                 Tolerance -> 10^-8, Shift -> currentShift}];
   (* form overlaps with last set of eigenvectors *)
   overlaps = Table[Abs[currentEvecs[[k]].#]& /@ eigensystem[[2]], {k, 3}]; 
   (* select α-continued eigenvectors *)
   maxOverlapPositions = Flatten[Position[#, Max[#]]& /@ overlaps];
   {currentEvals, currentEvecs} = #[[maxOverlapPositions]]& /@ eigensystem;
   (* store eigenvalues and overlaps *)
   evData[αExp] = {currentEvals, currentEvecs};
   overlapData[αExp] = {overlaps, maxOverlapPositions};
   (* make contour plots of eigenfunctions *)
   If[MemberQ[αExpValuesToShow, αExp],
    Show[GraphicsArray[eigenFunctionListContourPlot /@ currentEvecs],
        PlotLabel -> N[αExp]]],
   {αExp, 0, 3, 3/allFrames}]]];

animationFrames[5]

Here are some of the eigenfunctions for various α.

animationFrames[All]

The following two graphics show the α-dependence of the eigenvalues (left) and the  α-dependence of the overlaps for the three selected states. Overall the overlap integrals are nearly 1, but show clear dips around α-values where two eigenvalues become nearly degenerate.

With[{colors = {RGBColor[1, 0, 0], RGBColor[0, 1, 0], RGBColor[0, 0, 1]}},
 Show[GraphicsArray[
 {Graphics[Transpose[{colors,
             Line /@ Transpose[Table[{αExp, #}& /@ evData[αExp][[1]],
                                     {αExp, 0, 3, 3/60}]]}], Frame -> True],
  Graphics[Transpose[{colors,
            Line /@ Transpose[Table[{αExp, #}& /@ (Table[#1[[j]][[#2[[j]]]], 
            {j, 3}]& @@ overlapData[αExp]), {αExp, 0, 2, 3/60}]]}], 
           Frame -> True]}]]]

For the eigenfunctions of an equilateral triangle with Robin boundary conditions, see [✶]. For isospectral domains with mixed Dirichlet and Neumann boundary conditions, see [✶]. For general considerations of the mixed boundary value case, see [✶].