Warum gibt Mathematica für diese Gleichung falsche Eigenwerte an?

Nov 23 2020

Hier ist ein Eigenwertproblem in der Zylinderkoordinate: $$\mu(r)\frac{\partial}{\partial r} \left( \frac{1}{\mu(r)}\frac{1}{r}\frac{\partial (ru)}{\partial r} \right)=-p^2u$$Dabei ist p der erforderliche Eigenwert. Der Koeffizient ist$$\mu(r)=500, 0 \leq r \leq a_{1}\\ \mu(r)=1,a_{1}<r \leq a$$ mit $a_{1}=0.004,a=0.06$und die Randbedingung ist $$u(r=0)=0,\\ u(r=a)=0.$$ Mit dem Befehl "NDEigenvalues" und der Auswahl von "FiniteElement" habe ich folgende Codes geschrieben:

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

Dieser Code liefert die Antwort:

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

Die obigen Werte sind jedoch falsch. Tatsächlich kann dieses Problem mit den Bessel-Funktionen gelöst werden$J_{n}(x)$ und $Y_{n}(x)$. Mit diesem Analyseverfahren habe ich völlig unterschiedliche Eigenwerte gefunden:

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

Jetzt bin ich sicher, dass die mit der Analysemethode erhaltenen Werte korrekt sind (ich habe 1D FEM codiert, das die gleichen Ergebnisse wie die analytische liefert). Warum liefert der Befehl "NDEigenvalues" die falschen Ergebnisse?

ps: Einige Erklärungen für die Analysemethode. Das Problem wurde aus der Analyse des Magnetfeldes abgeleitet.$u(r)$ ist eine Komponente des Vektorpotentials.$\mu(r)$ist die relative Permeabilität. Daher sind Kontinuitäten an der Schnittstelle erforderlich. Wenn ich bezeichne$$u(r)=u_{1}(r), 0 \leq r \leq a_{1}\\ u(r)=u_{2}(r),a_{1}<r \leq a\\ \mu_{r}=500$$ Dann hätten wir haben sollen $$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}$$ Wenn ich dieses Problem mit der Analysemethode löse, kann ich zwei Ansätze für schreiben $u_{1}, u_{2}:$ $$u_{1}(r)=R_{1}(pa_{1})J_{1}(pr)\\ u_{2}(r)=J_{1}(pa_{1})R_{1}(pr)$$ Und die entsprechende Eigenwertgleichung ist $$\mu_{r}J_{1}(pa_{1})R_{0}(pa_{1})=J_{0}(pa_{1})R_{1}(pa_{1}) \quad (1)$$ wo $$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)$$Gl. (1) kann mit der Newton-Raphson-Methode gelöst werden, um die richtigen Eigenwerte zu erhalten.

Antworten

5 AlexTrounev Nov 24 2020 at 22:30

Dieses Problem in einem Fall des 3D-FEM-Vektorpotentials wird hier diskutiert . Wir können die Funktion approvon xzczd answer wie folgt verwenden

\[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}

Visualisierung

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

2 SPPearce Nov 29 2020 at 00:28

Ich habe ein Paket zum Lösen von 1D-Eigenwert-BVPs, das solche mit Schnittstellen enthält. Es konstruiert die "Evans-Funktion", eine analytische Funktion, die den Eigenwerten des ursprünglichen Systems entspricht, und reduziert das Problem darauf, Wurzeln einer glatten Funktion einer Variablen zu finden. Siehe meinen Github oder meine Antworten auf andere Fragen auf der Website.

Installieren Sie das Paket:

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

Wir müssen zuerst die resultierenden ODEs mit meiner Funktion ToMatrixSystem in eine Matrixform umwandeln:

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}

Dies hat noch einen nicht spezifizierten Wert $\epsilon$, der Grenzwert von $r \rightarrow 0$.

Für einen gegebenen Wert von $\epsilon$ und der Eigenwert $p$Wir können die Evans-Funktion bewerten. Zum Beispiel für$p=1$ und $\epsilon = 10^{-3}$::