Mathematicaがこの方程式に間違った固有値を与えるのはなぜですか？
円筒座標での固有値の問題は次のとおりです。 $$\mu(r)\frac{\partial}{\partial r} \left( \frac{1}{\mu(r)}\frac{1}{r}\frac{\partial (ru)}{\partial r} \right)=-p^2u$$ここで、pは必要な固有値です。係数は$$\mu(r)=500, 0 \leq r \leq a_{1}\\ \mu(r)=1,a_{1}<r \leq a$$ と $a_{1}=0.004,a=0.06$、および境界条件は $$u(r=0)=0,\\ u(r=a)=0.$$ 「NDEigenvalues」コマンドを使用して「FiniteElement」を選択すると、次のコードを記述しました。
μ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}
ただし、上記の値は正しくありません。実際、この問題はベッセル関数を使用して解決できます$J_{n}(x)$ そして $Y_{n}(x)$。この分析手順で、私はまったく異なる固有値を見つけました。
{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}
これで、分析メソッドによって取得された値が正しいことを確認できます（分析メソッドと同じ結果を提供する1D FEMをコーディングしました）。では、なぜ「NDEigenvalues」コマンドが間違った結果をもたらすのでしょうか。
ps：分析方法の説明。この問題は、磁場の分析から導き出されました。$u(r)$ ベクトルポテンシャルの成分です。$\mu(r)$は比透磁率です。したがって、インターフェースには連続性が必要です。私が示すなら$$u(r)=u_{1}(r), 0 \leq r \leq a_{1}\\ u(r)=u_{2}(r),a_{1}<r \leq a\\ \mu_{r}=500$$ その後、私たちは持っている必要があります $$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}$$ 分析法を使用してこの問題を解決するとき、私は2つの仮説を書くことができます $u_{1}, u_{2}:$ $$u_{1}(r)=R_{1}(pa_{1})J_{1}(pr)\\ u_{2}(r)=J_{1}(pa_{1})R_{1}(pr)$$ そして、対応する固有値方程式は次のとおりです。 $$\mu_{r}J_{1}(pa_{1})R_{0}(pa_{1})=J_{0}(pa_{1})R_{1}(pa_{1}) \quad (1)$$ どこ $$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)$$式 （1）はニュートン・ラフソン法で解き、正しい固有値を得ることができます。
回答
3D FEMベクトルポテンシャルの場合のこの問題については、ここで説明します。approxzczdanswerの関数は次のように使用できます
\[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}
視覚化
Table[Plot[fun[[i]], {r, 0, 1}, PlotLabel -> p[[i]]], {i, Length[p]}]
1D固有値BVPを解くためのパッケージがあります。これには、インターフェイスを持つものが含まれています。これは、元のシステムの固有値に対応する分析関数である「エバンス関数」を構築し、問題を1つの変数の滑らかな関数の根を見つけることに減らします。サイト上の私のgithubまたは他の質問への回答を参照してください。
パッケージをインストールします。
Needs["PacletManager`"]
PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]
まず、関数ToMatrixSystemを使用して、結果のODEを行列形式に変換する必要があります。
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}
これはまだ不特定の値を持っています $\epsilon$、の制限値 $r \rightarrow 0$。
の与えられた値に対して $\epsilon$ と固有値 $p$Evans関数を評価できます。たとえば、$p=1$ そして $\epsilon = 10^{-3}$：
Evans[1, sys /. ϵ -> 10^-3]
(* -1.53145*10^-6 *)
プロットは、この関数のいくつかのルーツがあることを示しています。
Plot[Evans[p, sys /. ϵ -> 10^-3], {p, 10, 200}]
そして、FindRoot特定の固有値を与えるために使用することができます：
FindRoot[Evans[p, sys /. ϵ -> 10^-3], {p, 10}]
(* {p -> 19.9443} *)
より高い精度のために、私たちは縮小することができます $\epsilon$ ゼロに向かって、そしてオプションをいじる：
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} *)