8
$\begingroup$

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}<r \leq a$$ 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]

This code provides the answer:

{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}<r \leq a\\ \mu_{r}=500$$ 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.

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

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 = Simplify`PWToUnitStep@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]}]

Figure 1

$\endgroup$
2
  • $\begingroup$ Thank you so much! It is very enlightening! $\endgroup$ – Karl Eichenhaus Nov 26 '20 at 16:18
  • $\begingroup$ @KarlEichenhaus You are welcome! $\endgroup$ – Alex Trounev Nov 26 '20 at 17:26
2
$\begingroup$

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}]

Plot of the Evans function, showing roots

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} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.