Pourquoi Mathematica donne des valeurs propres erronées pour cette équation?

Nov 23 2020

Voici un problème de valeur propre en coordonnée cylindrique: $$\mu(r)\frac{\partial}{\partial r} \left( \frac{1}{\mu(r)}\frac{1}{r}\frac{\partial (ru)}{\partial r} \right)=-p^2u$$où p est la valeur propre requise. Le coefficient est$$\mu(r)=500, 0 \leq r \leq a_{1}\\ \mu(r)=1,a_{1}<r \leq a$$ avec $a_{1}=0.004,a=0.06$, et la condition aux limites est $$u(r=0)=0,\\ u(r=a)=0.$$ En utilisant la commande "NDEigenvalues" et en choisissant "FiniteElement", j'ai écrit les codes suivants:

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

Ce code fournit la réponse:

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

Cependant, les valeurs ci-dessus sont incorrectes. En fait, ce problème peut être résolu en utilisant les fonctions de Bessel$J_{n}(x)$ et $Y_{n}(x)$. Avec cette procédure analytique, j'ai trouvé des valeurs propres totalement différentes:

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

Maintenant, je suis sûr que les valeurs obtenues par la méthode analytique sont correctes (j'ai codé 1D FEM qui fournit les mêmes résultats à l'analyse). Alors pourquoi la commande "NDEigenvalues" donne-t-elle de mauvais résultats?

ps: Quelques explications sur la méthode analytique. Le problème a été dérivé de l'analyse du champ magnétique.$u(r)$ est une composante du potentiel vectoriel.$\mu(r)$est la perméabilité relative. Par conséquent, des continuités sont requises sur l'interface. Si je dénote$$u(r)=u_{1}(r), 0 \leq r \leq a_{1}\\ u(r)=u_{2}(r),a_{1}<r \leq a\\ \mu_{r}=500$$ Alors nous devrions avoir $$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}$$ Lors de la résolution de ce problème en utilisant la méthode analytique, je peux écrire deux ansatzes pour $u_{1}, u_{2}:$ $$u_{1}(r)=R_{1}(pa_{1})J_{1}(pr)\\ u_{2}(r)=J_{1}(pa_{1})R_{1}(pr)$$ Et l'équation de valeur propre correspondante est $$\mu_{r}J_{1}(pa_{1})R_{0}(pa_{1})=J_{0}(pa_{1})R_{1}(pa_{1}) \quad (1)$$ où $$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) peut être résolu par la méthode de Newton-Raphson, pour obtenir les valeurs propres correctes.

Réponses

5 AlexTrounev Nov 24 2020 at 22:30

Ce problème dans un cas de potentiel vectoriel FEM 3D est discuté ici . Nous pouvons utiliser la fonction approde la réponse xzczd comme suit

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

2 SPPearce Nov 29 2020 at 00:28

J'ai un package pour résoudre les BVP à valeur propre 1D, qui inclut ceux avec des interfaces. Il construit la "Fonction Evans", une fonction analytique dont la correspondance correspond aux valeurs propres du système original, réduisant le problème à trouver les racines d'une fonction lisse d'une variable. Voir mon github ou mes réponses à d'autres questions sur le site.

Installez le package:

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

nous devons d'abord transformer les ODE résultants en une forme matricielle en utilisant ma fonction 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}

Cela a toujours une valeur non spécifiée $\epsilon$, la valeur limite de $r \rightarrow 0$.

Pour une valeur donnée de $\epsilon$ et la valeur propre $p$nous pouvons évaluer la fonction d'Evans. Par exemple, pour$p=1$ et $\epsilon = 10^{-3}$: