Dlaczego Mathematica podaje błędne wartości własne dla tego równania?

Nov 23 2020

Oto problem z wartością własną we współrzędnych cylindrycznych: $$\mu(r)\frac{\partial}{\partial r} \left( \frac{1}{\mu(r)}\frac{1}{r}\frac{\partial (ru)}{\partial r} \right)=-p^2u$$gdzie p jest wymaganą wartością własną. Współczynnik wynosi$$\mu(r)=500, 0 \leq r \leq a_{1}\\ \mu(r)=1,a_{1}<r \leq a$$ z $a_{1}=0.004,a=0.06$a warunek brzegowy to $$u(r=0)=0,\\ u(r=a)=0.$$ Używając polecenia „NDEigenvalues” i wybierając „FiniteElement”, napisałem następujące kody:

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

Ten kod dostarcza odpowiedzi:

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

Jednak powyższe wartości są nieprawidłowe. W rzeczywistości ten problem można rozwiązać za pomocą funkcji Bessela$J_{n}(x)$ i $Y_{n}(x)$. Dzięki tej procedurze analitycznej znalazłem zupełnie inne wartości własne:

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

Teraz jestem pewien, że wartości uzyskane metodą analityczną są poprawne (zakodowałem 1D MES, który daje takie same wyniki jak analityczny). Dlaczego więc polecenie „NDEigenvalues” daje błędne wyniki?

ps: Kilka wyjaśnień dotyczących metody analitycznej. Problem został wyprowadzony z analizy pola magnetycznego.$u(r)$ jest składową potencjału wektora.$\mu(r)$jest względną przepuszczalnością. W związku z tym w interfejsie wymagane są ciągłości. Jeśli oznaczę$$u(r)=u_{1}(r), 0 \leq r \leq a_{1}\\ u(r)=u_{2}(r),a_{1}<r \leq a\\ \mu_{r}=500$$ Wtedy powinniśmy $$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}$$ Rozwiązując ten problem metodą analityczną, mogę napisać dwie odpowiedzi $u_{1}, u_{2}:$ $$u_{1}(r)=R_{1}(pa_{1})J_{1}(pr)\\ u_{2}(r)=J_{1}(pa_{1})R_{1}(pr)$$ A odpowiadające mu równanie wartości własnej to $$\mu_{r}J_{1}(pa_{1})R_{0}(pa_{1})=J_{0}(pa_{1})R_{1}(pa_{1}) \quad (1)$$ gdzie $$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) można rozwiązać metodą Newtona-Raphsona, aby uzyskać prawidłowe wartości własne.

Odpowiedzi

5 AlexTrounev Nov 24 2020 at 22:30

Ten problem w przypadku potencjału wektorowego MES 3D jest omawiany tutaj . Możemy użyć funkcji approz odpowiedzi xzczd w następujący sposób

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

Wyobrażanie sobie

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

2 SPPearce Nov 29 2020 at 00:28

Mam pakiet do rozwiązywania 1D wartości własnych BVP, który obejmuje te z interfejsami. Konstruuje „funkcję Evansa”, funkcję analityczną, która odpowiada wartościom własnym oryginalnego układu, redukując problem do znalezienia pierwiastków gładkiej funkcji jednej zmiennej. Zobacz mój github lub moje odpowiedzi na inne pytania w witrynie.