Mathematica가이 방정식에 대해 잘못된 고유 값을 제공하는 이유는 무엇입니까?

Nov 23 2020

다음은 원통형 좌표의 고유 값 문제입니다. $$\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}

그러나 위의 값은 올바르지 않습니다. 실제로이 문제는 Bessel 함수를 사용하여 해결할 수 있습니다.$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"명령이 잘못된 결과를 제공하는 이유는 무엇입니까?

추신 : 분석 방법에 대한 몇 가지 설명. 문제는 자기장 분석에서 비롯되었습니다.$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}$$ 분석 방법을 사용하여이 문제를 풀 때, 나는 두 개의 ansatze를 쓸 수 있습니다. $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) 올바른 고유 값을 얻기 위해 Newton-Raphson 방법으로 풀 수 있습니다.

답변

5 AlexTrounev Nov 24 2020 at 22:30

3D FEM 벡터 전위의 경우이 문제는 여기서 논의 됩니다 . approxzczd 답변의 함수 를 다음과 같이 사용할 수 있습니다.

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

2 SPPearce Nov 29 2020 at 00:28

인터페이스가있는 것을 포함하는 1D 고유 값 BVP를 해결하기위한 패키지가 있습니다. 원래 시스템의 고유 값에 해당하는 분석 함수 인 "Evans 함수"를 구성하여 한 변수의 평활 함수의 근을 찾는 문제를 줄입니다. 내 github 또는 사이트의 다른 질문 에 대한 답변을 참조하십시오 .