Hcfyk

I am trying to compute the eigenfunctions of an oblate spheroid (a=75 cm and b=60 cm) using Mathematica's FEM package and Chris' answer from here. Specifically, I am looking for eigenfrequencies around 433, 893, 913 and 2400 MGHz. Is there any way I could narrow my search besides getting all eigenfrequencies initially and then looking for the desired outcome which is impractical?

Here is my code for the first 4 eigenmodes:

Needs["NDSolve`FEM`"];



helmholzSolve3D[g_, numEigenToCompute_Integer, 

opts : OptionsPattern] := 

Module[{u, x, y, z, t, pde, dirichletCondition, mesh, boundaryMesh, 

nr, state, femdata, initBCs, methodData, initCoeffs, vd, sd, 

discretePDE, discreteBCs, load, stiffness, damping, pos, nDiri, 

numEigen, res, eigenValues, eigenVectors, 

evIF},



(*Discretize the region*)



If[Head[g] === ImplicitRegion || Head[g] === ParametricRegion, 

mesh = ToElementMesh[DiscretizeRegion[g, opts], opts], 

mesh = ToElementMesh[DiscretizeGraphics[g, opts], opts]];

boundaryMesh = ToBoundaryMesh[mesh];



(*Set up the PDE and boundary condition*)



pde = D[u[t, x, y, z], t] - Laplacian[u[t, x, y, z], {x, y, z}] + 

u[t, x, y, z] == 0;

dirichletCondition = DirichletCondition[u[t, x, y, z] == 0, True];

(*Pre-process the equations to obtain the FiniteElementData in 

StateData*)nr = ToNumericalRegion[mesh];

{state} = 

NDSolve`ProcessEquations[{pde, dirichletCondition, 

u[0, x, y, z] == 0}, u, {t, 0, 1}, Element[{x, y, z}, nr]];

femdata = state["FiniteElementData"];

initBCs = femdata["BoundaryConditionData"];

methodData = femdata["FEMMethodData"];

initCoeffs = femdata["PDECoefficientData"];



(*Set up the solution*)vd = methodData["VariableData"];



sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];



(*Discretize the PDE and boundary conditions*)



discretePDE = DiscretizePDE[initCoeffs, methodData, sd];

discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];



(*Extract the relevant matrices and deploy the boundary conditions*)



load = discretePDE["LoadVector"];

stiffness = discretePDE["StiffnessMatrix"];

damping = discretePDE["DampingMatrix"];

DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];



(*Set the number of eigenvalues ignoring the Dirichlet positions*)



pos = discreteBCs["DirichletMatrix"]["NonzeroPositions"][[All, 2]];

nDiri = Length[pos];

numEigen = numEigenToCompute + nDiri;



(*Solve the eigensystem*)



res = Eigensystem[{stiffness, damping}, -numEigen];

res = Reverse /@ res;

eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];

eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];

evIF = ElementMeshInterpolation[{mesh}, #] & /@ eigenVectors;



(*Return the relevant information*)



{eigenValues, evIF, mesh}]



{ev, if, mesh} = 

helmholzSolve3D[Ellipsoid[{0, 0, 0}, {0.75, 0.6, 0.6}], 4, 

MaxCellMeasure -> 0.025]



Table[

DensityPlot[

if[[i]][x, y, 0.1], {x, -1, 1}, {y, -1, 1},

RegionFunction -> Function[{x, y}, x^2/0.75^2 + y^2/0.6^2 < 1],

PlotLabel -> ev[i] ,

ColorFunction -> Hue,

PlotLegends -> Automatic

],

{i, 1, 4}

]

Any suggestions?

You could use something like this:

{vals, funs} = 

 NDEigensystem[{-Laplacian[u[x, y, z], {x, y, z}] + u[x, y, z], 

   DirichletCondition[u[x, y, z] == 0, True]}, u, 

  Element[{x, y, z}, mesh], 4, 

  Method -> {"Eigensystem" -> {"FEAST", "Interval" -> {425, 500}}}]



{{427.961, 428.783, 430.026, 430.156},...}

You may try Eigensystem with

Method -> {"FEAST", "Interval" -> {a, b}}

to search eigenvalue pairs within an interval. See the documentation of Eigensystem, Section "Methods", Subsection "FEAST" for more details.