# Why Mathematica gives wrong eigenvalues for this equation?

Here is an eigenvalue problem in cylindrical coordinate: $$\mu(r)\frac{\partial}{\partial r} \left( \frac{1}{\mu(r)}\frac{1}{r}\frac{\partial (ru)}{\partial r} \right)=-p^2u$$ where p is the required eigenvalue. The coefficient is $$\mu(r)=500, 0 \leq r \leq a_{1}\\ \mu(r)=1,a_{1} with $$a_{1}=0.004,a=0.06$$,and the boundary condition is $$u(r=0)=0,\\ u(r=a)=0.$$ Using the "NDEigenvalues" command and choosing "FiniteElement", I've written following codes:

μr = 500; a1 = 4/10^3; a = 6/10^2;
μ = With[{μm = μr, μa = 1}, If[0 <= r <= a1, μm, μa]];
ℒ = μ*D[(1/μ)*(1/r)*D[r*u[r], r], r];
ℬ = DirichletCondition[u[r] == 0, True];
vals = NDEigenvalues[{ℒ, ℬ}, u[r], {r, 0, a}, 30,
Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.0001, "MaxBoundaryCellMeasure"-> 0.00001, "MeshOrder" -> 2}}}];
p = Sqrt[-vals]


{63.861766132883865, 116.92644447823088, 169.55780223711812, 222.06153226109987, 274.51050083985103, 326.93097516766255, 379.3347396704956,
431.7278681218963, 484.113808910877, 536.4946651790507, 588.8717924983509, 641.2461039100476, 693.6182368779678, 745.988649959372,
798.3576814523224, 850.7255863929587, 903.0925606857338, 955.4587573010893, 1007.8242974270114, 1060.1892783147352, 1112.5537789108064,
1164.9178639705115, 1217.2815871087598, 1269.6449930975, 1322.0081196163815, 1374.3709986038718, 1426.733657310317, 1479.0961191278266,
1531.458404249732, 1583.8205301993034}


However,the above values are incorrect. In fact, this problem can be solved using the Bessel functions $$J_{n}(x)$$ and $$Y_{n}(x)$$. With this analytical procedure, I've found totally different eigenvalues:

{19.750686053012217, 79.50553925115048, 136.9291955924841, 193.73804196226334, 250.2908871563726, 306.70770650924777, 363.04222591866534,
419.3226661586999, 475.56541618908665, 531.7806506165634, 587.9749498993451, 644.1526020560387, 700.3161917251147, 756.4665699161246,
812.6015250490414, 868.7082899215693, 924.6790897037489, 957.8509197090044, 981.4684330754833, 1037.3301171523472, 1093.4113326541358,
1149.5170337175198, 1205.62883441715, 1261.7420635874469, 1317.8550029034939, 1373.9668072980996, 1430.0768539865803, 1486.1843801285418,
1542.287997723794, 1598.3843930403937}


Now I'm sure that the values obtained by analytical method are correct (I've coded 1D FEM which provides the same results to the analytical one). So why does the "NDEigenvalues" command give the wrong results?

ps: Some explanations for the analytical method. The problem was derived from the analysis of the magnetic field. $$u(r)$$ is a component of the vector potential.$$\mu(r)$$ is the relative permeability. Hence, continuities are required on the interface. If I denote $$u(r)=u_{1}(r), 0 \leq r \leq a_{1}\\ u(r)=u_{2}(r),a_{1} Then we should have $$u_{1}(r)=0, r=0\\ u_{2}(r)=0, r=a\\ u_{1}(r)=u_{2}(r), r=a_{1}\\ \frac{1}{\mu_{r}}\frac{\partial}{\partial r}(ru_{1})=\frac{\partial}{\partial r}(ru_{2}),r=a_{1}$$ When solving this problem using the analytical method, I can write two ansatzes for $$u_{1}, u_{2}:$$ $$u_{1}(r)=R_{1}(pa_{1})J_{1}(pr)\\ u_{2}(r)=J_{1}(pa_{1})R_{1}(pr)$$ And the corresponding eigenvalue equation is $$\mu_{r}J_{1}(pa_{1})R_{0}(pa_{1})=J_{0}(pa_{1})R_{1}(pa_{1}) \quad (1)$$ where $$R_{1}(pr)=J_{1}(pr)Y_{1}(pa)-J_{1}(pa)Y_{1}(pr)\\ R_{0}(pr)=J_{0}(pr)Y_{1}(pa)-J_{1}(pa)Y_{0}(pr)$$ Eq. (1) can be solved by the Newton-Raphson method, to get the correct eigenvalues.

• Could you show your analytical solution code as well? Nov 23, 2020 at 11:43
• Please see the edited version. Nov 23, 2020 at 12:39
• It looks like you do not apply any boundary condition at $r=a$, while in numerical model there is DirichletCondition[u[r] == 0, True] Nov 23, 2020 at 15:34
• You can see that $u_{2}(a)=0$ is already satisfied, with the provided ansatz. And I've updated the description again. Thank you. Nov 24, 2020 at 9:23
• Ok! See my answer with solution of the problem. It is not exactly what you got, but data not so different. Nov 24, 2020 at 19:16

This problem in a case of 3D FEM vector potential is discussed here. We can use function appro from xzczd answer as follows

\[Mu]r = 500; a1 = 4/10^3; a = 6/10^2; d = a1/a;
\[Mu] = With[{\[Mu]m = \[Mu]r, \[Mu]a = 1},
If[0 <= r <= d, \[Mu]m, \[Mu]a]]; appro =
With[{k = 2 10^5}, ArcTan[k #]/Pi + 1/2 &];
mu = SimplifyPWToUnitStep@PiecewiseExpand@If[r <= d, \[Mu]r, 1] /.
UnitStep -> appro;
\[ScriptCapitalL] = mu D[1/mu (1/r)*D[r*u[r], r], r]/a^2;
\[ScriptCapitalB] = DirichletCondition[u[r] == 0, True];
{vals, fun} =
NDEigensystem[{\[ScriptCapitalL], \[ScriptCapitalB]},
u[r], {r, 0, 1}, 10,
Method -> {"PDEDiscretization" -> {"FiniteElement", {"MeshOptions" \
-> {"MaxCellMeasure" -> 0.00001}}}}];
p = Sqrt[-vals]

Out[]= {19.9785, 79.8404, 137.385, 194.307, 250.965, 307.482, 363.911, 420.282, 476.611, 532.91}


Visualisation

Table[Plot[fun[[i]], {r, 0, 1}, PlotLabel -> p[[i]]], {i, Length[p]}]


• Thank you so much! It is very enlightening! Nov 26, 2020 at 16:18
• @KarlEichenhaus You are welcome! Nov 26, 2020 at 17:26

I have a package for solving 1D eigenvalue BVPs, that includes those with interfaces. It constructs the "Evans Function", an analytic function whose correspond to the eigenvalues of the original system, reducing the problem to finding roots of a smooth function of one variable. See my github or my answers to other questions on the site.

Install the package:

Needs["PacletManager"]
PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]


we first need to turn the resulting ODEs into a matrix form using my function ToMatrixSystem:

sys = ToMatrixSystem[{D[1/r  D[r u1[r], r], r] + p^2 u1[r] == 0,
D[1/r  D[r u2[r], r], r] + p^2 u2[r] == 0},
{u1[ϵ] == 0, u2[a] == 0, u1[a1] == u2[a1],
1/μr (D[r u1[r], r] /. r -> a1) == (D[r u2[r], r] /. r -> a1) },
{u1, u2}, {r, ϵ, a1, a}, p] /. {μr -> 500, a1 -> 4/10^3, a -> 6/10^2}


This still has an unspecified value $$\epsilon$$, the limiting value of $$r \rightarrow 0$$.

For a given value of $$\epsilon$$ and the eigenvalue $$p$$ we can evaluate the Evans function. For instance, for $$p=1$$ and $$\epsilon = 10^{-3}$$:

Evans[1, sys /. ϵ -> 10^-3]
(* -1.53145*10^-6 *)


A plot shows there are some roots of this function:

Plot[Evans[p, sys /. ϵ -> 10^-3], {p, 10, 200}]


And then FindRoot can be used to give specific eigenvalues:

FindRoot[Evans[p, sys /. ϵ -> 10^-3], {p, 10}]
(* {p -> 19.9443} *)


For higher precision, we can shrink $$\epsilon$$ towards zero, and fiddle with the options:

p /. FindRoot[Evans[p, sys /. ϵ -> 10^-10, NormalizationConstants -> {0, 1},
WorkingPrecision -> 50], {p, #}, WorkingPrecision -> 50] & /@ {10, 100, 150, 200} // Quiet

(* {19.7506836087553767185196899913,
79.5055392302968147610410441291,
136.929195538974955894770829013,
193.738041724568292657607041215,
250.290886522212012980557959916} *)