क्यों गणितज्ञ इस समीकरण के लिए गलत eigenvalues ​​देता है?

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$$जहां पी को आवश्यक स्वदेशी है। गुणांक है$$\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}

अब मुझे यकीन है कि विश्लेषणात्मक विधि द्वारा प्राप्त मूल्य सही हैं (मैंने 1 डी एफईएम कोडित किया है जो विश्लेषणात्मक एक ही परिणाम प्रदान करता है)। तो "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}$$ विश्लेषणात्मक पद्धति का उपयोग करके इस समस्या को हल करते समय, मैं दो ansatzes लिख सकता हूं $u_{1}, u_{2}:$ $$u_{1}(r)=R_{1}(pa_{1})J_{1}(pr)\\ u_{2}(r)=J_{1}(pa_{1})R_{1}(pr)$$ और इसी eigenvalue समीकरण है $$\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)$$Eq। (1) न्यूटन-राफसन विधि द्वारा हल किया जा सकता है, ताकि सही स्वदेशी प्राप्त हो सके।

जवाब

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} 

VISUALIZATION

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

2 SPPearce Nov 29 2020 at 00:28

मेरे पास 1D eigenvalue BVPs को हल करने के लिए एक पैकेज है, जिसमें इंटरफेस वाले शामिल हैं। यह "इवांस फंक्शन" का निर्माण करता है, एक विश्लेषणात्मक फ़ंक्शन जिसका मूल सिस्टम के आइजनवेल्स के अनुरूप है, एक चर के एक चिकनी फ़ंक्शन की जड़ों को खोजने के लिए समस्या को कम करता है। साइट पर अन्य प्रश्नों के लिए मेरे गिथब या मेरे उत्तर देखें ।

पैकेज स्थापित करें:

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

हमें सबसे पहले अपने फ़ंक्शन ToMatrixSystem का उपयोग करके परिणामस्वरूप ODEs को मैट्रिक्स रूप में बदलना होगा:

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$हम इवांस फ़ंक्शन का मूल्यांकन कर सकते हैं। उदाहरण के लिए, के लिए$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} *)