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.

Numerics Volume, Chapter 1

Consider a square Koch curve of order 2. Solve the Helmholtz equation -Δ u(x,y)=(ω^2) u(x,y) in this Koch curve with Dirichlet boundary conditions. Estimate accuracy of the eigenvalues through comparison with state that is locally everywhere a product state of sine functions. Visualize some of the eigenfunctions as 3D plots and as contour plots. Visualize the values of the gradients of some of the lowest and some of the higher eigenfunctions.
Calculate the first 1000 lowest eigenvalues and compare the number of eigenfunctions n(ω) as function of ω with the asymptotic Weyl–Berry law [✶], [✶], [✶], [✶], [✶], [✶], [✶]

n(ω)~A/(4π) (ω^2)-{({{c (ω^d), for fractal boundaries}, {a /(4π)ω, for smooth boundaries.}})

The solve the eigenvalue problem of the Helmholtz equation, the piecewise horizontal and vertical boundary segments suggest a finite difference discretization of the interior of the Koch curve. To achieve a solution of sufficiently high precision, we will use higher order finite difference approximations.

Off[General::spell]; Off[General::spell1];

We start by calculating the 1D finite difference weights. For the actual 2D Koch curve problem, we use a tensor product form of these differences.

The kth element of FDWeightLists[ν][o] gives the list of length o of the weights of the function values to approximate the νth derivative at the kth point from the left through the function values at o points. It is obtained by forming an interpolating polynomial through the (symbolic) function values, differentiating this polynomial, and extracting the coefficients of the function values [✶], [✶].

FDWeightLists[ν_][o_] :=
Module[{ipo, poly}, 
       (* interpolating polynomial *) 
       ipo[x_] = InterpolatingPolynomial[Table[{x0 + k h, f[k]}, 
                                               {k, 0, o - 1}], x];
       (* loop over all o grid points *)
       Table[(* value of derivative at j-th grid point *)
           poly = Expand[h^ν Together[D[ipo[x], {x, ν}] /. x -> x0 + j h]];
           (* coefficients of function values *)
           Coefficient[poly, Table[f[k], {k, 0, o - 1}]], {j, 0, o - 1}]]

Here are two examples: the first and second derivative of a stencil of length three.

FDWeightLists[1][3]

FDWeightLists[2][3]

The function hOrder[ν, fdws] calculates the order of the approximation of the derivative (f^((ν)))(x) for the finite difference weights fdws at the point x and step size h.

hOrder[ν_, fdws_, {x_, h_}] := 
Normal[(* keep only highest order in h *) 
       Series[#, {h, 0, Position[Rest[#[[3]]], _?(# =!= 0&), 
            {1}, 1, Heads -> False][[1, 1]]}]]& /@ 
Series[(* expand in series about h == 0 *)
(* form expression approximating the derivative *)
MapIndexed[Function[{w, p}, w.Table[f[x + k h], 
               {k, 1 - p[[1]], -p[[1]] + Length[w]}]/h^ν], fdws],
       {h, 0, Length[fdws[[1]]]}]

Here are the approximation properties of the second derivative using three and five points. We see that the derivatives at the center points are approximated to an higher order.

hOrder[2, FDWeightLists[2][3], {x, h}]

hOrder[2, FDWeightLists[2][5], {x, h}]

We remark that the finite difference weights can be obtained from one of the built-in supporting functions from the NDSolve` context, namely through the function NDSolve`FiniteDifferenceDerivative.

FDWeightListsBI[ν_][o_] :=
With[{grid = Range[o]}, Coefficient[#, f /@ grid]& /@ 
     NDSolve`FiniteDifferenceDerivative[Derivative[ν], grid, 
                        f /@ grid, DifferenceOrder -> o - ν]]

(* quick check that the two functions give identical weights *)
Table[FDWeightLists[ν][o] == FDWeightListsBI[ν][o], 
      {ν, 1, 6}, {o, ν + 1, 6}] // Flatten // Union

The following graphic shows the central finite difference approximation weights for the second derivative. The graphic shows substantial changes for the first few values of α (we take a stencil of size 2α+1). As a result, we expect substantially different results for the eigenvalues for small α. Later, we will use α=3, which will give results accurate to less than 0.1% over a large ω-range.

Show[Graphics[Reverse @ Table[{Hue[(α - 1)/10], 
Line[MapIndexed[{#2[[1]] - α - 1, #1}&,
                FDWeightListsBI[2][2 α + 1][[α + 1]]]]},
                {α, 7}]], PlotRange -> All, Frame -> True]

For later convenience, we define one more finite difference weight-related function. Given a point on a 1D grid, the function stencilWeightList gives the list of weights approximating the vth derivative by using leftPoints to the left and rightPoints to the right (and the center point itself).

stencilWeightList[ν_][leftPoints_, rightPoints_] :=
        stencilWeightList[ν][leftPoints, rightPoints] = 
        FDWeightListsBI[ν][leftPoints + 1 + rightPoints][[leftPoints + 1]]

Next, we construct the Koch drum. A Koch drum is obtained by recursively applying the following rule edgeRule to the line segments of the last iteration drum. The 0th order drum is just a square.

edgeRule[{p_, q_}] := 
Module[{d = (q - p)/4, n = -Reverse[q - p]{-1, 1}/4},
       Partition[FoldList[Plus, p, {d, n, d, -n, -n, d, n, d}], 2, 1]]
       
startSquare = 
     Partition[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}, {-1, -1}}, 2, 1];

The following graphic shows the zeroth order drum in red, the first order in green, and the second order in violet. In the following, we will calculate the eigenvibrations of the second order Koch drum. The underlying square grid that will be used to partition the interior of the drum in the following is shown too.

(* list of the first few the Koch drum iterations *)
iterationOrder = 2;
nl = NestList[Flatten[edgeRule /@ #, 1]&, startSquare, iterationOrder];

(* edge length *)
λ = 2^-(2 iterationOrder - 1);

Module[{L = 2 - 2^-iterationOrder}, 
Show[Graphics[{
 {(* background mesh *)
  Thickness[0.001], GrayLevel[0.8], {#, Map[Reverse, #, {3}]}&[
   Table[Line[{{-L, z}, {L, z}}], {z, -L, L, λ}]]},
   MapIndexed[{Hue[0.78 (#2[[1]] - 1)/
               If[iterationOrder == 0, 1, iterationOrder]],   
  Thickness[0.01/#2[[1]]], Line /@ #1}&, nl]}], 
  PlotRange -> All, AspectRatio -> Automatic, Axes -> True]]

These are the line segments forming the boundary of the Koch drum. To test if a line segment is part of the boundary, we implement the function outerEdgeQ.

drumEdges = Sort /@ nl[[-1]];
(outerEdgeQ[#] = True)& /@ drumEdges;

translateSquare[mp_, edge_, square_]  :=
With[{t = 2 ((Plus @@ edge/2) - mp)}, Map[(t + #)&, square, {-2}]]

Until now, we only have the boundary of the Koch drum. While for a boundary element method, this would be enough to calculate the eigenvalues, for a finite difference approximation, we need explicit points in the interior. The interior of the Koch drum is naturally subdivided into squares of edge lengths equal to the length of a boundary segments. We obtain all squares by starting in the center and recursively select all neighboring squares until we hit the boundary. The list allInnerSquares contains the inner squares in the form Square[midpoint, squareEdges].

(* form squares that are neighbors of a given square sq *)
neighborSquares[sq:Square[mp_, edges_]] :=
Module[{sortedEdges = Sort /@ edges, growEdges, newSquares}, 
       (* the edges of the square that are not boundaries of the drum *)
       growEdges = Complement[sortedEdges, Select[sortedEdges, outerEdgeQ]];
       newSquares = translateSquare[mp, #, sq]& /@ growEdges]

allInnerSquares = 
Module[{usedMpQ, currentSquares, squaresBag, newSquareLayer},
If[(* zeroth order is special *)
 iterationOrder === 0, {Square[{0, 0}, startSquare]},
 (* the four squares around the center *)
 currentSquares = Square[Plus @@ Flatten[#, 1]/8, #]& /@  
  (Partition[Append[#, First[#]], 2, 1]& /@ 
     ({{{ 0,  0}, { λ,  0}, { λ,  λ}, { 0,  λ}}, 
       {{-λ,  0}, { 0,  0}, { 0,  λ}, {-λ,  λ}},
       {{ 0, -λ}, { λ, -λ}, { λ,  0}, { 0,  0}}, 
       {{-λ, -λ}, { 0, -λ}, { 0,  0}, {-λ,  0}}}));
(* currently already used squares *)
 squaresBag = currentSquares;
 (usedMpQ[#[[1]]] = True)& /@ currentSquares;
(* form outside neighbor square layers as long as possible *)
While[newSquareLayer = First /@ Split[Sort[Flatten[neighborSquares /@
                                                   currentSquares]],
                                      #1[[1]] === #2[[1]]&];
      currentSquares = Complement[newSquareLayer, 
                                  Select[newSquareLayer, usedMpQ[#[[1]]]&]];
      currentSquares =!= {},
      (* update *)
      squaresBag = {squaresBag, currentSquares};
     (usedMpQ[#[[1]]] = True)& /@ currentSquares];
(* the squares forming the interior of the Koch drum *)
allInnerSquares = Flatten[squaresBag]]];

Here are all the 256 interior squares.

Show[Graphics[{(* the interior squares, colored randomly *)
               {Hue[Random[]], #}& /@ allInnerSquares /. 
                Square[mp_, edges_] :> Polygon[First /@ edges],
               (* the drum edges *)
               {GrayLevel[0], Thickness[0.005], Line /@ drumEdges}},
      AspectRatio -> Automatic]]

We generate lists of all vertices, the outer (boundary) vertices, and the strictly inner vertices.

allVertices = Union[Flatten[Last /@ allInnerSquares, 2]];
outerVertices = Union[Flatten[drumEdges, 1]];
innerVertices = Complement[allVertices, outerVertices];

To obtain a good resolution of the solution in the thin bays of the Koch drum, we will now subdivide each square into (rL^2) smaller squares. We use rL=8 in the following. L denotes the shortest distance between two neighboring discretization points.

(* order of refinement within each square *) rL = 8;
(* distance between two adjacent vertices *) L = λ/rL;

(* subdivide a square into rL×rL smaller squares *)
refineSquare[Square[mp_, edges_], rL_]  :=
Module[{p1, p2, p3, p4, d1, d2},
       {p1, p2, p3, p4} = First /@ edges;
       d1 = p2 - p1; d2 = p4 - p1;
       Table[Square[{p1 + j/rL d1 + k/rL d2, p1 + (j + 1)/rL d1 + k/rL d2,
                     p1 + (j + 1)/rL d1 + (k + 1)/rL d2, 
                     p1 + j/rL d1 + (k + 1)/rL d2}],
             {j, 0, rL - 1}, {k, 0, rL - 1}] // Flatten];

We now create lists of inner vertices innerRefinedVertices (of the small squares) and vertices boundaryVertices (of the small) squares that lie on the boundary of the drum.

allRefinedInnerSquares = Flatten[refineSquare[#, rL]& /@ allInnerSquares];
allRefinedVertices = Union[Flatten[First /@ allRefinedInnerSquares, 1]];

boundaryVertices = 
Union[Flatten[Table[#1 + j/rL (#2 - #1), {j, 0, rL}]& @@@ drumEdges, 1]];

innerRefinedVertices = Complement[allRefinedVertices, boundaryVertices];

So we have 15361 inner vertices. At these inner vertices, we will have to calculate the values of the eigenvibrations. (Due to the Dirichlet boundary conditions we are using all boundary vertex values are identically zero.)

Length /@ {allRefinedVertices, boundaryVertices, innerRefinedVertices}

Here is a sketch of the Koch drum with the inner vertices in red and the boundary vertices in black.

Show[Graphics[{{PointSize[0.004], GrayLevel[0], Point /@ boundaryVertices},
               {PointSize[0.002], Hue[0], Point /@ innerRefinedVertices}}],
     PlotRange -> All, AspectRatio -> Automatic]

Now the values of the eigenvibrations at the inner vertices have to be determined. This means that we have to find eigenvalues and eigenvectors of a d×d matrix, where for the Dirichlet boundary conditions d is the number of the interior vertices.

dimD = Length[innerRefinedVertices];

We define some functions that number the vertices, return their coordinates, and test if a vertex is an inner vertex.

(* convert between vertex coordinates and vertex number *)
Do[vertexDNumber[innerRefinedVertices[[k]]] = k;
   vertexDCoordinates[k] = innerRefinedVertices[[k]];
   isInnerDVertex[innerRefinedVertices[[k]]] = True, {k, dimD}];
   
(* is a vertex an inner vertex? *)
isInnerDVertexQ[vertex_] := isInnerDVertexQ[vertex] = 
                            TrueQ[isInnerDVertex[vertex]]

We define similar functions for the boundary vertices.

Do[boundaryVertexNumber[boundaryVertices[[k]]] = k;
   boundaryVertexCoordinates[k] = boundaryVertices[[k]];
   isBoundaryVertex[boundaryVertices[[k]]] = True, 
   {k, Length[boundaryVertices]}];
   
(* is a vertex a boundary vertex? *)
isBoundaryVertexQ[vertex_] := isBoundaryVertexQ[vertex] = 
                              TrueQ[isBoundaryVertex[vertex]]

For all inner vertices, we now define a neighborhood and approximate the Laplacian through the weighted average of the function values at neighboring vertices. The auxiliary function addVertexWeights adds the contributions of the list of finite difference weights vws for each vertex.

addVertexWeights[vws_] := 
If[Length[#] === 1, First[#], #[[1, 1]] -> (Plus @@ (Last /@ #))]& /@ 
                              Split[Sort[vws], #1[[1]] === #2[[1]]&]

For a given stencil of odd length 2α+1, the function makeNonVanishingDMatrixElements returns a list of the nonvanishing elements of the discretization of the negative of the 2D Laplace operator. We incorporate the Dirichlet boundary conditions through an antisymmetric continuation of the function across the boundary. By using an odd stencil length, we make sure that is always possible.

makeNonVanishingDMatrixElements[α_Integer] := addVertexWeights @
Module[{centerVertex, jB, swl = -stencilWeightList[2][α, α]},
Table[centerVertex = innerRefinedVertices[[k]]; 
       {{k, k} -> 2 swl[[α + 1]],
        Function[dir, jB =.;
        (* loop over neighbors *)
        Table[nextVertex = centerVertex + j L dir;
              (* neighbor is inner vertex, boundary vertex, 
                 or outer vertex *)
              Which[isInnerDVertexQ[nextVertex], 
                    {k, vertexDNumber[nextVertex]} -> swl[[α + 1 + j]], 
                    isBoundaryVertexQ[nextVertex], jB = j; Sequence @@ {},
                    True, reflectedVertex = (* reflect *)
                             centerVertex + jB L dir - (j - jB) L dir;
                    {k, vertexDNumber[reflectedVertex]} -> 
                                        -swl[[α + 1 + j]]], 
                    {j, 1, α}]] /@ (* horizontal and vertical neighbors *)
                    {{+1, 0}, {-1, 0}, {0, 1}, {0, -1}}}, 
                    (* all inner vertices *) {k, dimD}] // Flatten];

Here is the resulting sparse matrix for stencil length seven. MD7 is of dimensions 15361×15361 and has 187429 nonvanishing elements. By construction, the matrix MD7 is symmetric.

MD7 = SparseArray[makeNonVanishingDMatrixElements[3]]

(* check that MD7 is a symmetric matrix *)
Max[MD7 - Transpose[MD7]]

We visualize the nonzero elements of the resulting matrix. It is nicely banded.

Show[Graphics[{GrayLevel[0], Rectangle[#1 - 1/2, #1 + 1/2]}& @@@ 
                                    makeNonVanishingDMatrixElements[3]],
     PlotRange -> All, Frame -> True, AspectRatio -> Automatic]

To find the eigenvalues and eigenvectors of the matrix MD7, we will use the "Arnoldi" method. We obtain the extremal eigenvalues with the "BothEnds" criteria. We multiply by the square of the discretization length L to obtain the eigenvalues for the original drum.

Eigenvalues[N[MD7], 2, Method -> {"Arnoldi", Criteria -> "BothEnds",  
                                  BasisSize -> 30}]/L^2

Next, we calculate the eight lowest eigenvalues. The function drumEigensystem calculates some selected eigenvalues of the matrix MD7. We see the unique lowest eigenvalue and that the next two eigenvalues are degenerate. Because of the fourfold rotational symmetry of the Koch curve, such a degeneracy is expected.

drumEigensystem[n_, opts___] := MapAt[1/L^2 #&, Reverse /@ 
    Eigensystem[N[MD7], n, Method -> {"Arnoldi", opts, Tolerance -> 10^-10,
                BasisSize -> Abs[20 n], MaxIterations -> 10^4}], 1]

{evalsD7Low, evecsD7Low} = drumEigensystem[-8];
ListPlot[evalsD7Low, PlotRange -> All, PlotStyle -> {PointSize[0.03]},
         AxesOrigin -> {0, 0}]

We also calculate the eigenvalues near the exactly known eigenvalues of the state, that is in each subsquare a product state of sine functions [✶]. The lowest of these functions has the eigenvalue Λ=2(((π/λ))^2)≈1263.31. The following list of eigenvalues clearly suggests that the third element of evecsD7Sin1 is this state. Its eigenvalue deviates by only 0.0006% from the exact value.

Λ = 2 (Pi/λ)^2;
{evalsD7Sin1, evecsD7Sin1} = drumEigensystem[6, Shift -> Λ L^2];
(evalsD7Sin1/Λ - 1)

The next state that is a product of sine functions has eigenvalue 4Λ. The eigenvalue of the fourth state of the list deviates by about 0.04% from the exact value. This indicates that the states up to this energy are clearly precisely enough calculated for any visualization purposes we are primarily interested in.

{evalsD7Sin2, evecsD7Sin2} = drumEigensystem[6, Shift -> 4 Λ L^2];
(evalsD7Sin2/(4 Λ) - 1)

Now that we have the eigenvectors, we will visualize them. We start with 3D plots. To construct them, we just use the small squares allRefinedInnerSquares from above and add as the z-value the eigenvector components to their vertices.

(* plot a general surface with height values over the drum vertices *)
drum3DPlot[vertexValueFunction_, opts___] :=
Graphics3D[{EdgeForm[], (* add z-component from eigenvector evec *)
            Map[vertexValueFunction, 
            (* the small squares *)             
            Polygon /@ (List @@ (First /@ allRefinedInnerSquares)), {-2}],
                {GrayLevel[0], Thickness[0.006], (* drum boundary *)
                 Line /@ Map[Append[#, 0]&, drumEdges, {-2}]}},
          opts, PlotRange -> All, BoxRatios -> {1, 1, 0.4}];

(* plot Dirichlet eigenfunctions *)
eigenDVibration3DPlot[evec_, opts___] :=
With[{ev = evec/Sign[Plus @@ evec]},
     drum3DPlot[Function[v, Append[v, If[isInnerDVertexQ[v], 
                            ev[[vertexDNumber[v]]], 0]]], opts]]

Here are the eigenvibrations of the nine lowest states.

Show[GraphicsArray[#]]& /@
     Partition[eigenDVibration3DPlot /@ evecsD7Low, 2];

The 3D plots of the higher eigenvalues show generally complicated patterns of the minima and maxima. The fourth of the above states though, is very regular. It is the lowest state that is globally a product of sine functions.

Show[GraphicsArray[#]]& /@
     Partition[eigenDVibration3DPlot /@ evecsD7Sin1, 2]

Here are two higher eigenfunctions whose oscillating parts are largely concentrated in the four bay areas [✶].

Show[GraphicsArray[eigenDVibration3DPlot[
       drumEigensystem[1, Shift -> # L^2][[2, 1]]]& /@ {176.642, 525.582}]]

As an alternative to the last 3D plots and as a better visualization method for the highly oscillating higher states, we will now implement a contour plot for the Koch drum. The built-in contour plot will, of course, not work directly. Two possibilities to construct a contour plot for the semiregular shape of the Koch drum come to mind: a) Extending the Koch curve to a square array by adding vertices with value zero, and b) making independent contour plots for all of the small squares and displaying them together. We will follow the method b) in the following. (This method is more flexible and can be easily generalized to say, the triangular Koch curve) The function eigenVibrationContourPlot generates a contour graphic of the eigenvalue evec using contours contours.

(* vertex numbers inside the refined squares  *)
refinedSquareVertexArray[
 Square[mp_, {{p1_, _}, {p2_, _}, {p3_, _}, {p4_, _}}], rL_]  :=
Module[{d1 = p2 - p1, d2 = p4 - p1, vertex},
       {Table[vertex = p1 + j/rL d1 + k/rL d2;
              If[isInnerDVertexQ[vertex], vertexDNumber[vertex], 0],
              {k, 0, rL}, {j, 0, rL}],
        {Min[#], Max[#]}& /@ Transpose[{p1, p2, p3, p4}]}];
        
refinedSquareData = refinedSquareVertexArray[#, rL]& /@ allInnerSquares;

(* contour values for contours *) 
contourValues[evec_, contours_] := 
Module[{eD = Length[evec]}, (* homogeneous contour spacing *)
       #[[Round[eD/contours/2]]]& /@ 
           Partition[Sort[evec], Round[eD/contours]]]

(* contour plot inside each of the squares of the Koch drum *)
smallSquareContourPlot[evec_, {vertexArray_, meshRange_}, cvs_] :=
Module[{lcpData},
       lcpData = Map[If[# =!= 0, evec[[#]], 0.]&, vertexArray, {2}];
       Graphics[ListContourPlot[lcpData, 
       (* set options so that the single plots fit together *)
        Contours -> cvs, MeshRange -> meshRange, PlotRange -> {0, 1},
        ColorFunctionScaling -> False, ContourLines -> False, 
        ContourSmoothing -> True, ColorFunction -> GrayLevel]]]

eigenVibrationContourPlot[evec_, contours_, opts___] :=
Module[{ev = evec/Sign[Plus @@ evec], min, max, cs},
       min  = Min[ev]; max = Max[ev];
       (* scale eigenvector values to [0, 1] *)
       scaledEv = (ev - min)/(max - min);
       cvs = contourValues[scaledEv, contours];
       Block[{$DisplayFunction = Identity},
        (* make contour plots in all small squares *)
        cs = smallSquareContourPlot[scaledEv, #, cvs]& /@ refinedSquareData;
        Show[{cs, (* the Koch drum boundary *)
              Graphics[{{Hue[0], Thickness[0.008], 
                             Line /@ drumEdges}}]},
             PlotRange -> All, Frame -> False]] /. 
                               Identity -> $DisplayFunction]

Here are contour plots for the nine lowest states.

Show[GraphicsArray[#]]& /@
     Partition[eigenVibrationContourPlot[#, 30]& /@ evecsD7Low,  2]

Show[GraphicsArray[#]]& /@
     Partition[eigenVibrationContourPlot[#, 30]& /@ evecsD7Sin1, 2]

Show[GraphicsArray[#]]& /@
     Partition[eigenVibrationContourPlot[#, 30]& /@ evecsD7Sin2, 2]

Now let us look at the Weyl–Berry law. A direct call to Eigenvalues[N[MD7], -1000, …] is not realistic for most computers due to the large memory requirements (about 10 GB for a basis size of 10 times the number of eigenvalues). So, instead we will move through the spectrum and calculate 24 eigenvalues at once, align the newly calculated eigenvalues with the last batch and extract the new higher eigenvalues. Because some eigenvalues have finite multiplicity, we will align four eigenvalues. Using this method we will get the first 1000 states in about an hour with the moderate memory usage of less than 300 MB. Because the Arnoldi method gives not only the eigenvalues, but the eigenvectors too, we also form the cumulative sum of the squares of eigenvectors evecSum. And we keep a truncated version of evecSum, called evecSumSmall that sums only states for which ω≳2π/λ where λ is the length of a straight boundary segment of our Koch curve.

(* center of the next set of eigenvalues *)
nextEvalCenter[l_] := (Plus @@ Take[l, -4]/4) + 
                      (* avoid eigenvalue  *) 1/2 (l[[5]] - l[[4]])

evalNumberMax = 1000;

(* start data *)
{currentEvals, currentEvecs} = {evalsD7Low, evecsD7Low};
evalList = currentEvals; evecSum = Plus @@ (currentEvecs^2); 


While[
 (* get additional, higher eigenvalues *)
 evDim = 24;
 While[(* the next higher eigenvalues and corresponding vectors *)
 {nextEvals, nextEvecs} = 
      drumEigensystem[evDim, Shift -> L^2 nextEvalCenter[evalList]];
 (* align with last set of eigenvalues *)
 δs = #.#& /@ ((Take[evalList, {-4, -1}] - #)& /@  
                                Partition[nextEvals, 4, 1]);
 minPos = Position[δs, _?((# < 10^-6)&)];                        
 (* genuinely new eigenvalues *)
 newPosis = If[minPos =!= {{Length[δs]}} && minPos =!= {}, 
               {minPos[[1, 1]] + 4, Length[nextEvals]}, {}];
 newPosis === {}, (* in case of spectral gaps use additive shift *)  
                  evDim = Round[3/2 evDim]];
 (* update list of eigenvalues and sum of eigenvectors squared *)
 evalList = Join[evalList, Take[nextEvals, newPosis]];
 evecSum = evecSum + Plus @@ (Take[nextEvecs, newPosis]^2);
 If[Max[nextEvals] < (Pi/λ)^2, evecSumSmall = evecSum];
 (* use evalNumberMax lowest states *)
 Length[evalList] < evalNumberMax, Null]

Length[evalList]

Here are the cumulative sums of the squares of the eigenfunctions for the first 135 and all 1000 eigenfunctions. These castles over the Koch curve show a flat Thomas–Fermi interior density and pronounced spikes at the boundaries.

Show[GraphicsArray[
 Block[{DisplayFunction = Identity},
       eigenDVibration3DPlot[#, Boxed -> False]& /@ 
                               {evecSumSmall, evecSum}]]]

For the visualization of the Weyl–Berry law, we subtract the leading term A/(4π) (ω^2) to obtain a less steep curve.

{drumBoundaryLength, drumArea} = 
  {Plus @@ (Sqrt[#.#]& /@ (Subtract @@@ drumEdges)),
   Plus @@ (Sqrt[#.#]^2& /@ (Subtract @@@ (#[[2, 1]]& /@ allInnerSquares)))}

The following graphic shows that the c (ω^d) term is a good approximation for ω≲0.9 2π/λ and the a /(4π)ω term for larger ω [✶]. The fitted exponent agrees well with fractal dimension 1.5 of the Koch drum.

(* subtract volume term *)
ΔIDOS = MapIndexed[{#1, -(#2[[1]] - 1/Pi #1^2)}&, Sqrt[evalList]];
(* fit exponent *)
FindFit[Select[ΔIDOS, #[[1]] < 22&], a x^b, {{a, 0.6}, {b, 1.5}}, x]

Show[Graphics[{
{Hue[0], Thickness[0.002], Line[Table[{ω, 0.57 ω^(3/2)}, {ω, 0, 40}]]},
{Hue[0.8], Thickness[0.002], Line[Table[{ω, 32/(4Pi) ω}, {ω, 0, 60}]]},
{GrayLevel[0], PointSize[0.006], Point /@ ΔIDOS}}], 
     Frame -> True, PlotRange -> All, FrameLabel -> {"ω", "δIDOS"}]

For a boundary with more corners, we should take into account contributions 1/(128 (π^2)) for π/2 corners and -5/(1152(π^2)) for 3/2π corners [✶].

We continue with the visualization of the absolute values of the gradients for the lowest eigenstates. We use again the finite difference approximations derived above, this time for the first derivative. We now use noncentered finite difference formulas.

grad1DValue[vertex_, dir_, maxStencilWidth_, evec_] :=
Module[{jL, jR, j, α, sData, swl, vcs}, 
(* special rule for boundary vertices *)
If[isBoundaryVertexQ[vertex],
   Which[isBoundaryVertexQ[vertex + L dir] || 
         isBoundaryVertexQ[vertex - L dir], {jL, jR} = {0, 0},
         isInnerDVertexQ[vertex + L dir], {jL, jR} = {0, maxStencilWidth}, 
         isInnerDVertexQ[vertex - L dir], {jL, jR} = {maxStencilWidth, 0}],
(* go to the left and the right of an inner vertex *)
{jL, jR} = (j = 1; While[nextVertex = vertex + # j L dir;
                         j <= maxStencilWidth + 1 && 
                         isInnerDVertexQ[nextVertex], j++]; 
            j - 1)& /@ {-1, +1};     
 α = (maxStencilWidth - 1)/2; 
 (* use maxStencilPoints *)
 Which[jL > α && jR > α, {jL, jR} = {α, α},
       jL + jR > maxStencilWidth,
       If[jR >= jL, 
          jL = Min[α, jL]; jR = maxStencilWidth - jL,
          jR = Min[α, jR]; jL = maxStencilWidth - jR]]]; 
 (* the weight vector *)
 swl = stencilWeightList[1][jL, jR];
 (* approximate value of the 1D derivative *)
 vcs = Table[vertexDNumber[vertex + j L dir], {j, -jL, jR}] /.
                                           _vertexDNumber :> 0;
 (* deal special with boundary vertices *)
 Which[{jL, jR} === {0, 0}, 0, 
       vcs[[+1]] === 0, Rest[swl].evec[[Rest[vcs]]],
       vcs[[-1]] === 0, Most[swl].evec[[Most[vcs]]],
       True, swl.evec[[vcs]]]]

gradValue[vertex_, maxStencilWidth_, evec_] := Sqrt[#.#]& @
{grad1DValue[vertex, {1, 0}, maxStencilWidth, evec],
 grad1DValue[vertex, {0, 1}, maxStencilWidth, evec]}

Given a list of the absolute values of the gradient vector at the inner vertices of the Koch drum, the function eigenVibrationGrad3DPlot makes a 3D plot of the values over the Koch drum domain.

eigenVibrationGrad3DPlot[gradInnerVector_, gradBoundaryVector_, opts___] :=
     drum3DPlot[Function[v, Append[v, If[isInnerDVertexQ[v], 
                     gradInnerVector[[vertexDNumber[v]]], 
                     gradBoundaryVector[[boundaryVertexNumber[v]]]]]]]

Here are the magnitudes of the gradients for the four lowest states. The gradients have singularities due to the 270° corners of the Koch curve at some points [✶], [✶], [✶], [✶] (for better discretizations near corners, see [✶], [✶]).

(* turn of message for low order finite differences *)
Off[NDSolve`FiniteDifferenceDerivative::aord];
Off[NDSolve`FiniteDifferenceDerivative::nonopt];

Show[GraphicsArray[#]]& /@ 
Partition[eigenVibrationGrad3DPlot[
 Table[gradValue[vertexDCoordinates[k], 7, evecsD7Low[[#]]], {k, dimD}],
 Table[gradValue[boundaryVertexCoordinates[k], 7, evecsD7Low[[#]]], 
       {k, Length[boundaryVertices]}]]& /@ {1, 2, 4, 5}, 2]

We end with a small excursion into the world of eigenfunctions with Neumann boundary conditions [✶], [✶]. Now we have to take into account all boundary vertices too. This yields 17409 vertices.

dimN = Length[allRefinedVertices];

Similar to the Dirichlet case, we define various functions for numbering and checking the type of a vertex.

Do[vertexNNumber[allRefinedVertices[[k]]] = k;
   vertexNCoordinates[k] = allRefinedVertices[[k]];
   isInnerNVertex[allRefinedVertices[[k]]] = True, {k, dimN}];
   
isInnerNVertexQ[vertex_] := isInnerNVertexQ[vertex] = 
                            TrueQ[isInnerNVertex[vertex]]

This time, we also differentiate between horizontal and vertical edges.

horizontalEdges = Select[drumEdges, #[[1, 2]] == #[[2, 2]]&];
verticalEdges   = Select[drumEdges, #[[1, 1]] == #[[2, 1]]&];

{horizontalBoundaryVertices, verticalBoundaryVertices} = 
Union[Flatten[Table[#1 + j/rL (#2 - #1), {j, 0, rL}]& @@@ #, 1]]& /@ 
{horizontalEdges, verticalEdges};

Do[isHorizontalBoundaryVertex[horizontalBoundaryVertices[[k]]] = True, 
   {k, Length[horizontalBoundaryVertices]}];
   
Do[isVerticalBoundaryVertex[verticalBoundaryVertices[[k]]] = True, 
   {k, Length[verticalBoundaryVertices]}];   
   
isHorizontalBoundaryVertexQ[vertex_] := isHorizontalBoundaryVertexQ[vertex] = 
TrueQ[isHorizontalBoundaryVertex[vertex]]

isVerticalBoundaryVertexQ[vertex_] := isVerticalBoundaryVertexQ[vertex] = 
TrueQ[isVerticalBoundaryVertex[vertex]]

isBoundaryVertexQ[dir_, vertex_] := isBoundaryVertexQ[dir, vertex] = 
Which[Abs[dir] === {1, 0}, isVerticalBoundaryVertexQ[vertex],
      Abs[dir] === {0, 1}, isHorizontalBoundaryVertexQ[vertex]]

For the construction of the matrix to be diagonalized, we proceed largely as in the Dirichlet case. But this time, we continue the function symmetrically across the boundary vertices to ensure Neumann boundary conditions. We again use a stencil of length 7.

nonVanishingNMatrixElements = addVertexWeights @
Module[{maxStencilSize = 7, α, swl, centerVertex, jB},
α = (maxStencilSize - 1)/2;
swl = -stencilWeightList[2][α, α];
(tab = Table[centerVertex = allRefinedVertices[[k]]; 
       {{k, k} -> 2 swl[[α + 1]],
        Function[dir, jB =.;
        If[isBoundaryVertexQ[dir, centerVertex], jB = 0];
        Table[nextVertex = centerVertex + j L dir; 
              Which[isInnerNVertexQ[nextVertex], 
                    If[isBoundaryVertexQ[dir, nextVertex], jB = j];
                    {k, vertexNNumber[nextVertex]} -> swl[[α + 1 + j]],
                    True, (* continue symmetrically *)
                    reflectedVertex = centerVertex + jB L dir - (j - jB) L dir;
                    {k, vertexNNumber[reflectedVertex]} -> 
                                        +swl[[α + 1 + j]]], 
                    {j, 1, α}]] /@ (* horizontal and vertical neighbors *)
                    {{+1, 0}, {-1, 0}, {0, 1}, {0, -1}}}, 
                    (* all inner vertices *) {k, dimN}]) // Flatten];

We form the sparse matrix of dimension 17409×17409 with 214005 nonvanishing elements.

MN7 = SparseArray[nonVanishingNMatrixElements]

Unfortunately, because of the symmetric continuation, the resulting matrix is not symmetric. (Using a discretization scheme that has discretization points shifted by L/2, so that the boundary would lie in the middle between two points would result in a symmetric matrix [✶].)

Max[MN7 - Transpose[MN7]]

The function drumNEigensystem calculates the eigenvalues and eigenfunctions of MN7. This time we use a higher tolerance because the nonsymmetric nature of the matrix will result in numerically complex eigenvalues and eigenfunctions.

drumNEigensystem[n_, opts___] := MapAt[1/L^2 #&, Reverse /@ 
    Eigensystem[N[MN7], n, Method -> {"Arnoldi", opts, Tolerance -> 10^-12,
                BasisSize -> 20 n, MaxIterations -> 10^4}], 1]

Here are a few of the states. The lowest eigenvalue is basically zero, reflecting the u(x,y)=constant solution. As expected, some eigenvalues have small imaginary parts.

{evalsN7Low, evecsN7Low} = drumNEigensystem[12, Shift -> 0]; 
evalsN7Low

We visualize the eigenfunctions using the routine eigenNVibration3DPlot. This time, we do not show the drum outline, because the functions do not vanish on the boundary.

(* plot Neumann eigenfunctions *)
eigenNVibration3DPlot[evec_, opts___] :=
With[{ev = Re[evec]}, DeleteCases[#, _Line, Infinity]& @
     drum3DPlot[Function[v, Append[v, ev[[vertexNNumber[v]]]]], opts,
                Boxed -> False, PlotRange -> 1.01 Max[Abs[ev]] {-1, 1}]]

Here are the first few eigenfunctions.

Show[GraphicsArray[#]]& /@ Partition[
             eigenNVibration3DPlot /@ evecsN7Low[[{1, 2, 4, 5, 6, 12}]], 3]

And here are three slightly higher states that are prominently localized in the four bays.

Show[GraphicsArray[eigenNVibration3DPlot /@
              drumNEigensystem[3, Shift -> 98.8646 L^2][[2]]]]

We end here and leave it to the reader to explore more of the world of eigenvalues and eigenfunctions of 2D domains with fractal boundaries. One could now continue and investigate the distribution of the eigenvalue spacings, and the distribution of the eigenfunction values [✶], [✶], [✶], [✶]. One could calculate and analyze the eigensystems of triangular Koch drums [✶] (which are slightly less interesting, because the bay areas are less pronounced and so the eigenfunctions can not reach them easily) and other shapes (for instance the areas bounded by some of the space-filling Peano and Hilbert curves). (For experimental results on fractal Koch drums, see [✶], [✶]; for classical paths in square Koch drums, see [✶].)