डिफरेंशियल क्वैडचर विधि नॉनलाइन बीसी के साथ 4 डी ऑर्डर पीडीई पर विफल हो जाती है क्योंकि ग्रिड सघन हो जाती है
मैं निम्नलिखित प्रारंभिक-सीमा मूल्य समस्या को हल करने के लिए अंतर चतुर्भुज विधि (DQM) का उपयोग कर रहा हूं:
$$\frac{\partial^2 w}{\partial t^2} + \mu \frac{\partial w}{\partial t} + C_0 \frac{\partial^4 w}{\partial x^4}=0$$ $$w(0,t)=X_b(t)$$ $$\frac{\partial w}{\partial x}\bigg|_{x=0}=\frac{\partial^2 w}{\partial x^2}\bigg|_{x=1}=\frac{\partial^3 w}{\partial x^3}\bigg|_{x=1}=0$$
यहाँ $X_b(t)$( Xb[t]नीचे दिए गए कोड में) सिस्टम में इनपुट है। हार्मोनिक समारोह के लिए
$$X_b(t)=2G_1 \cos(2\Omega t)$$
इनपुट के रूप में, NDSolveपूरी तरह से काम करता है। अन्य इनपुट के लिए भी, सिमुलेशन ठीक से चलता है। लेकिन इनपुट के लिए
$$X_b(t)=2G \cos(2 \Omega t) (w(1,t)+\alpha \; w(1,t)^3)$$
उच्च आवृत्ति दोलनों में वृद्धि होती है और सिमुलेशन परिणाम अधिक से अधिक गलत हो जाते हैं, क्योंकि ग्रिड बिंदुओं की संख्या ( Npनीचे दिए गए कोड में) बढ़ जाती है। यदि tmaxबड़ा है या Np > 10, तो सिमुलेशन चेतावनी के साथ विफल हो जाता है
NDSolve :: ndsz: एकवचन या कठोर प्रणाली संदिग्ध।
मैंने इस समस्या को एक अन्य विधि से हल किया है, कोई उच्च आवृत्ति दोलनों मौजूद नहीं है।
निम्नलिखित मेरा परीक्षण है। पीडीई को Np - 1DQM के साथ ODEs में दिया गया है ।
Np = 10; G1 = 0.05; Ω = 1; μ = 0.05; tmax = 10; a = 30;
ii = Range[1, Np]; X = 0.5 (1 - Cos[(ii - 1)/(Np - 1) π]);
Xi[xi_] := Cases[X, Except[xi]];
M[xi_] := Product[xi - Xi[xi][[l]], {l, 1, Np - 1}];
C1 = C3 = ConstantArray[0, {Np, Np, 4}];
Table[If[i != j, C1[[i, j, 1]] = M[X[[i]]]/((X[[i]] - X[[j]]) M[X[[j]]])];,
{i, 1, Np}, {j, 1, Np}];
Table[C1[[i, i, 1]] = -Total[Cases[C1[[i, All, 1]], Except[C1[[i, i, 1]]]]];,
{i, 1, Np}];
C3[[All, All, 1]] = C1[[All, All, 1]];
C3[[1, All, 1]] = 0;
C3[[All, All, 2]] = C1[[All, All, 1]].C3[[All, All, 1]];
C3[[Np, All, 2]] = 0;
C3[[All, All, 3]] = C1[[All, All, 1]].C3[[All, All, 2]];
C3[[Np, All, 3]] = 0;
C3[[All, All, 4]] = C1[[All, All, 1]].C3[[All, All, 3]];
C3r4 = N[C3[[All, All, 4]]];
C0 = 1.8751^-4;
K1M = C0 C3r4[[2 ;; Np, 1 ;; Np]];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, Np}];
α = -0.001;
इनपुट Xb[t]को स्तंभ वेक्टर के माध्यम से समीकरणों की प्रणाली में प्रतिस्थापित किया जाता है YV[t]।
Xb[t] = 2 G1 Cos[2 Ω t] (Subscript[x, Np][t] + α Subscript[x, Np][t]^3);
YV[t] = Flatten[{Xb[t], Y1[t]}];
eqns = Thread[D[Y1[t], t, t] + μ D[Y1[t], t] + K1M.YV[t] == 0];
Lg = 1; bt = 1.8751/Lg; ξ = X[[2 ;; Np]];
y0 = -0.5 a (((Cos[bt*ξ] - Cosh[bt*ξ])-0.734*(Sin[bt*ξ] - Sinh[bt*ξ])));
X0 = -y0; X0d = 0 y0;
s = NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax},
Method -> {"StiffnessSwitching", Method->{"ExplicitRungeKutta", Automatic}},
MaxSteps -> Infinity,
AccuracyGoal -> 11, PrecisionGoal -> 11]; //AbsoluteTiming
plot1 = Plot[Evaluate[Subscript[x, Np][t] /. First@s], {t, 0, tmax},
PlotRange -> All]
11.1 के लिए संस्करण के साथ प्राप्त परिणाम Np = 6और 8नीचे दिए गए हैं। के लिए Np = 8, आउटपुट में उतार-चढ़ाव बढ़ता है।
के लिए Np = 12, यह चेतावनी देता है
NDSolve :: ndsz: t == 7.12986041201692828 पर, चरण आकार प्रभावी रूप से शून्य है; एकवचन या कठोर प्रणाली का संदेह।
मैंने NDSolveकठोर प्रणाली से निपटने के लिए विभिन्न विकल्पों का उपयोग किया है, लेकिन फिर भी यह काम नहीं कर रहा है।
मैंने सोचा कि, अगर मैं NDSolveछोटे अंतराल में उपयोग करता हूं तो यह सफल हो सकता है। इसलिए मैंने कोड बनाया जिसमें प्रारंभिक स्थितियों (ith पुनरावृत्ति के लिए) की गणना पिछले पुनरावृत्ति (i - 1) से आउटपुट के आधार पर की जाती है। लेकिन कई NDSolveसिमुलेशन भी काम नहीं किया।
मैंने विभिन्न प्रारंभिक स्थितियों में भी कोशिश की है, लेकिन चेतावनी बनी हुई है। किसी ने मुझे इस मुद्दे को हल करने में मदद कर सकते हैं? धन्यवाद।
जवाब
यदि DQM को सही तरीके से कार्यान्वित किया जाता है, तो यह विधि की एक आवश्यक सीमा हो सकती है। मुझे DQM के बारे में कुछ भी नहीं पता था, लेकिन इस पेपर को स्कैन करने पर , मुझे लगता है कि यह विधि समान है Pseudospectral। दरअसल, एक त्वरित परीक्षण से पता चलता है कि, DQM में व्युत्पन्न 1 ऑर्डर व्युत्पन्न का भार गुणांक मैट्रिक्स बिल्कुल वैसा ही है जैसा कि pseudospectral विधि में 1 क्रम व्युत्पन्न के विभेदन मैट्रिक्स:
NDSolve`FiniteDifferenceDerivative[1, X, DifferenceOrder -> "Pseudospectral"][
"DifferentiationMatrix"] == C1[[All, All, 1]]
(* True *)
हालांकि मैं इस समय एक ठोस उदाहरण नहीं दे सकता, लेकिन मैंने देखा कि Pseudospectralजब कुछ मामलों में स्थानिक ग्रिड बिंदु बढ़ जाते हैं तो अस्थिर हो सकता है। यदि आपकी समस्या इस प्रकार की है तो परीक्षण करें। विशेष bc की वजह से, हम NDSolveसमस्या को हल करने के लिए सीधे उपयोग नहीं कर सकते हैं , तो चलिए सिस्टम को अलग करते हैं$x$अपने आप से दिशा। मैं pdetoodeकार्य के लिए उपयोग करेंगे ।
G1 = 0.05; Ω = 1; μ = 0.05; tmax = 10; a = 30; α = -0.001;
Lg = 1; bt = 1.8751/Lg; C0 = 1.8751^-4;
icrhs = 0.5 a ((Cos[bt x] - Cosh[bt x]) - 0.734 (Sin[bt x] - Sinh[bt x]));
With[{w = w[x, t]},
eqn = D[w, t, t] + μ D[w, t] + C0 D[w, {x, 4}] == 0;
bc = {(w /. x -> 0) == (2 G1 Cos[2 Ω t] (w + α w^3) /. x -> 1),
D[w, x] == 0 /. x -> 0,
{D[w, x, x] == 0, D[w, {x, 3}] == 0} /. x -> 1};
ic = {w == icrhs, D[w, t] == 0} /. t -> 0];
Off[NDSolve`FiniteDifferenceDerivative::ordred]
domain = {0, 1};
CGLGrid[{xl_, xr_}, n_] := 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)]);
del = #[[3 ;; -3]] &;
help[points_] := (grid = Array[# &, points, domain];
grid = N@CGLGrid[domain, points];
difforder = points;
(*Definition of pdetoode isn't included in this post,
please find it in the link above.*)
ptoofunc = pdetoode[w[x, t], t, grid, difforder];
ode = ptoofunc[eqn] // del;
odeic = del /@ ptoofunc[ic];
odebc = {ptoofunc@bc[[1, 1]] == ptoofunc@bc[[1, -1]], ptoofunc@bc[[2 ;;]]};
NDSolveValue[{ode, odeic, odebc}, w@grid[[-1]], {t, 0, tmax}, MaxSteps -> ∞])
wRlst = Table[help[points], {points, 6, 10}];
ListLinePlot /@ wRlst
ListLinePlot[help[11], PlotRange -> {-30, 30}]
NDSolveValue :: ndsz: at t == 2.356194489774355`, स्टेप साइज़ प्रभावी रूप से शून्य है; एकवचन या कठोर प्रणाली का संदेह।
जैसा कि हम देख सकते हैं, जब ग्रिड बिंदुओं की संख्या 10 से अधिक नहीं होती है, तो समाधान स्थिर प्रतीत होता है और काफी तेजी से परिवर्तित होता है, लेकिन एक बार pointsवृद्धि होने के बाद 11, समाधान ओपी के कार्यान्वयन के व्यवहार के समान जंगली हो जाता है।
तो, कैसे करें चक्कर? कम क्रम के अंतर सूत्र का उपयोग करने के लिए विवेकशील होना प्रभावी होता है:
points = 50; grid = Array[# &, points, domain];
difforder = 2;
(* Definition of pdetoode isn't included in this post,
please find it in the link above. *)
ptoofunc = pdetoode[w[x, t], t, grid, difforder];
ode = ptoofunc[eqn] // del;
odeic = del /@ ptoofunc[ic];
odebc = {ptoofunc@bc[[1, 1]] == ptoofunc@bc[[1, -1]], ptoofunc@bc[[2 ;;]]};
tst = NDSolveValue[{ode, odeic, odebc}, w@grid[[-1]], {t, 0, tmax}, MaxSteps -> ∞];
ListLinePlot[tst, PlotRange -> {-30, 30}]
जैसा कि ऊपर दिखाया गया है, समाधान के साथ स्थिर रहता है points = 50; difforder = 2।
difforder = 4 अगर आप की तरह भी इस्तेमाल किया जा सकता है।
परिशिष्ट: ओपी की विधि का पुन: कार्यान्वयन
इस उत्तर की शुरुआत में पेपर के साथ ओपी के कोड को एक साथ देखने के बाद, मुझे लगता है कि मैंने समझ लिया है कि ओपी ने क्या लागू किया है। निम्नलिखित उसी विधि के लिए मेरा कार्यान्वयन है, जो थोड़ा अधिक पठनीय हो सकता है:
G1 = 0.05; Ω = 1; μ = 0.05; a = 30; C0 = 1/1.8751^4; Lg = 1; bt = 1.8751/Lg; α = -0.001;
tmax = 10;
domain = {0, 1};
CGLGrid[{xl_, xr_}, n_] := 1/2 (xl + xr + (xl - xr) Cos[(π Range[0, n - 1])/(n - 1)]);
points = 9;
grid = N@CGLGrid[domain, points];
d1 = NDSolve`FiniteDifferenceDerivative[1, grid, DifferenceOrder -> "Pseudospectral"][
"DifferentiationMatrix"];
(* Alternative method for defining d1: *)
(*
Π[i_] := Times@@Delete[grid[[i]] - grid, {i}];
d1 = Table[If[i == k,0.,Π@i/((grid[[i]] - grid[[k]]) Π@k)], {i, points}, {k, points}];
Table[d1[[i, i]] = -Total@d1[[i]], {i, points}];
*)
d1withbc = Module[{d = d1}, d[[1, All]] = 0.; d];
d2 = d1 . d1withbc;
d2withbc = Module[{d = d2}, d[[-1, All]] = 0.; d];
d3 = d1 . d2withbc;
d3withbc = Module[{d = d3}, d[[-1, All]] = 0.; d];
d4 = d1 . d3withbc;
W[t_] = Rest@Array[w[#][t] &, points];
w0 = 2 G1 Cos[2 Ω t] (w[points][t] + α w[points][t]^3);
eqns = W''[t] + μ W'[t] + C0 Rest[d4] . Flatten@{w0, W[t]} == 0 // Thread;
ξ = Rest@grid;
X0 = 0.5 a ((Cos[bt ξ] - Cosh[bt ξ]) - 0.734 (Sin[bt ξ] - Sinh[bt ξ]));
X0d = 0 ξ;
sol = NDSolve[{eqns, W[0] == X0, W'[0] == X0d}, W[t], {t, 0, tmax},
MaxSteps -> ∞][[1]]; // AbsoluteTiming
Plot[w[points][t] /. sol // Evaluate, {t, 0, 10}, PlotRange -> All]
कुछ और स्पष्टीकरण: इस विधि में $\frac{\partial^4}{\partial x^4}$ को पुनरावर्ती तरीके से अलग कर दिया गया है $\frac{\partial}{\partial x}$पहले ( C1[[All, All, 1]]ओपी के कोड में, d1मेरे कोड में) और विवेकाधीन है$\frac{\partial^4}{\partial x^4}$का उपयोग करके गणना की जाती है Dot। अभी भी संदिग्ध लग रहा है? ठीक है, चलो मान्य करें:
f[x_] = (x - 1/2)^6 + Sin[x];
ListPlot[{grid, #}\[Transpose] & /@ {f'[grid], d1.f[grid]}, PlotMarkers -> {"o", "x"},
PlotLegends -> {"exact f'", "approximated f'"}]
ListPlot[{grid, #}\[Transpose] & /@ {f''[grid], d1.d1.f[grid]},
PlotMarkers -> {"o", "x"}, PlotLegends -> {"exact f''", "approximated f''"}]
ListPlot[{grid, #}\[Transpose] & /@ {f'''[grid], d1.d1.d1.f[grid]},
PlotMarkers -> {"o", "x"}, PlotLegends -> {"exact f'''", "approximated f'''"}]
जबसे $\frac{\partial}{\partial x}$, $\frac{\partial^2}{\partial x^2}$ तथा $\frac{\partial^3}{\partial x^3}$ सभी विधि में मध्यवर्ती के रूप में प्रकट हुए हैं, ओपी की समस्या के बीएसी को सीधे मैट्रिक्स को संशोधित करके लगाया जा सकता है, उदाहरण के लिए:
ListLinePlot[{grid, d1withbc . f[grid]}\[Transpose]]
जैसा कि ऊपर दिखाया गया है, $\frac{\partial f}{\partial x}\Big|_{x=0}=0$ लगाया गया है।
चूंकि यह कोड ब्रैकट बीम के लिए DQM लागू कर रहा है, इसलिए हमें इस कोड को ग्रिड Npपरिवर्तन की संख्या के साथ स्थिर बनाने के लिए सही सीमा स्थिति रखने की आवश्यकता है । यह केवल छोटा संशोधन है लेकिन यह किसी भी Npउदाहरण के लिए काम करता है
Np = 20; G1 = 0.05; Ω = 1; μ = 0.05; tmax = 10; a = 30;
ii = Range[1, Np]; X = 0.5 (1 - Cos[(ii - 1)/(Np - 1) π]);
Xi[xi_] := Cases[X, Except[xi]];
M[xi_] := Product[xi - Xi[xi][[l]], {l, 1, Np - 1}]; C1 =
C3 = ConstantArray[0, {Np, Np, 4}];
Table[If[i != j,
C1[[i, j, 1]] = M[X[[i]]]/((X[[i]] - X[[j]]) M[X[[j]]])];, {i, 1,
Np}, {j, 1, Np}];
Table[C1[[i, i,
1]] = -Total[Cases[C1[[i, All, 1]], Except[C1[[i, i, 1]]]]];, {i, 1, Np}];
C3[[All, All, 1]] = C1[[All, All, 1]];
C3[[1, All, 1]] = 0 C1[[1, All, 1]];
C3[[All, All, 2]] = C1[[All, All, 1]].C3[[All, All, 1]];
C3[[Np, All, 2]] = 0 C1[[Np, All, 2]];
C3[[All, All, 3]] = C1[[All, All, 1]].C3[[All, All, 2]];
C3[[Np, All, 3]] = 0 C1[[Np, All, 2]];
C1[[All, All, 2]] = C1[[All, All, 1]].C1[[All, All, 1]];
C3[[All, All, 4]] = C1[[All, All, 1]].C3[[All, All, 3]];
C3r4 = N[C3[[All, All, 4]]]; C11 = C3r4[[1, 1]]; C0 = 1.8751^-4;
K1M = C0 C3r4[[2 ;; Np, 1 ;; Np]]; K1V = C0 C3r4[[2 ;; Np, 1]];
Y1[t_] := Table[Subscript[x, i][t], {i, 2, Np}];
a2 = Flatten[ConstantArray[1, {Np - 1, 1}]]; α = -0.001;
Xb[t] = 2 G1 Cos[2 Ω t] (Subscript[x, Np][
t] + α Subscript[x, Np][t]^3) Table[KroneckerDelta[Np, i], {i, 2, Np}];
YV[t] = Flatten[{0, Y1[t]}];
eqns = Thread[D[Y1[t], t, t] + μ D[Y1[t], t] + K1M.YV[t] == Xb[t]];
Lg = 1; bt = 1.8751/Lg; ξ = X[[2 ;; Np]];
y0 = -0.5 a (((Cos[bt*ξ] - Cosh[bt*ξ]) -
0.734*(Sin[bt*ξ] - Sinh[bt*ξ]))); X0 = -y0; X0d = 0 y0;
s = NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d},
Y1[t], {t, 0, tmax}]; // AbsoluteTiming
plot1 = Plot[Evaluate[Subscript[x, Np][t] /. First@s], {t, 0, tmax},
PlotRange -> All, PlotLabel -> Row[{"Np = ", Np}]]
इस दृष्टिकोण के साथ हम इस पर विचार करने के लिए है Xb[t]के रूप में बाहरी बल मनमाने ढंग से बात करने के लिए लागू किया जाता nextके रूप में
Xb[t] = 2 G1 Cos[
2 Ω t] (Subscript[x, Np][
t] + α Subscript[x, Np][t]^3) Table[
KroneckerDelta[next, i], {i, 2, Np}];
मामले में next=Npहम ऊपर कोड है। मूल कोड अस्थिर समाधान का उत्पादन करने का मुख्य कारण है K1M, Hamed Sadeghian, Ghaz Rezazadeh और Peter M. Osterberg द्वारा MEMS स्विच के Pull-In Phenomena के सामान्यीकृत विभेदक चतुष्कोण विधि के अनुप्रयोग से ली गई मैट्रिक्स परिभाषा। हम इसे इस प्रकार पुनर्परिभाषित कर सकते हैं
Np = 10; G1 = .05; \[CapitalOmega] = 1; \[Mu] = 0.05; tmax = 10; a = \
30;
ii = Range[1, Np]; X = 0.5 (1 \[Minus] Cos[(ii - 1)/(Np - 1) \[Pi]]);
Xi[xi_] := Cases[X, Except[xi]];
M[xi_] := Product[xi - Xi[xi][[l]], {l, 1, Np - 1}]; C1 =
ConstantArray[0, {Np, Np}];
Table[If[i != j,
C1[[i, j]] = M[X[[i]]]/((X[[i]] - X[[j]]) M[X[[j]]])];, {i,
Np}, {j, Np}];
Table[C1[[i,
i]] = -Total[Cases[C1[[i, All]], Except[C1[[i, i]]]]];, {i, Np}];
W1 = C1; W1[[1, All]] = 0.;
W2 = C1.W1; W2[[Np, All]] = 0.;
W3 = C1.W2; W3[[Np, All]] = 0.;
W4 = C1.W3; W4[[1, All]] = 0.;
C0 = 1.8751^-4; K1M = C0 W4;
Y1[t_] := Table[Subscript[x, i][t], {i, 1, Np}]; Y4 = K1M.Y1[t];
\[Alpha] = -0.001; Xb =
2 G1 Cos[2 \[CapitalOmega] t] (Subscript[x, Np][
t] + \[Alpha] Subscript[x, Np][t]^3) Table[
KroneckerDelta[1, i], {i, 1, Np}];
eqns = Thread[D[Y1[t], t, t] + \[Mu] D[Y1[t], t] + Y4 == Xb];
Lg = 1; bt = 1.8751/Lg; \[Xi] = X;
y0 = -0.5 a (((Cos[bt*\[Xi]] - Cosh[bt*\[Xi]]) -
0.734*(Sin[bt*\[Xi]] - Sinh[bt*\[Xi]]))); X0 = -y0; X0d = 0 y0;
s = NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax}];
हम इस समाधान की तुलना Xb=0लाल बिंदुओं पर points=10(ठोस बिंदु) xzczd कोड द्वारा उत्पन्न समाधान से कर सकते हैं
अब यदि हम ऊपर दिए गए कोड की तरह पहले बिंदु पर भी Np=30लागू होते हैं Xb, तो हमें निम्न प्रकार से प्रत्येक ग्रिड बिंदु के लिए चित्र मिलता है
Table[Plot[Evaluate[Subscript[x, i][t] /. First@s], {t, 0, tmax},
PlotRange -> All, PlotLabel -> Row[{"i = ", i}]], {i, 1, Np}]
यह बाहरी बल के लिए बहुत आम प्रतिक्रिया है। इस मैट्रिक्स का उपयोग करके K1M = C0 W4हम Xbकार्यान्वयन के मुख्य विचार को महसूस कर सकते हैं$x_1(t)$ निम्नलिखित नुसार
Np = 12; G1 = .05; \[CapitalOmega] = 1; \[Mu] = 0.05; tmax = 20; a = \
30;
ii = Range[1, Np]; X = 0.5 (1 \[Minus] Cos[(ii - 1)/(Np - 1) \[Pi]]);
Xi[xi_] := Cases[X, Except[xi]];
M[xi_] := Product[xi - Xi[xi][[l]], {l, 1, Np - 1}]; C1 =
ConstantArray[0, {Np, Np}];
Table[If[i != j,
C1[[i, j]] = M[X[[i]]]/((X[[i]] - X[[j]]) M[X[[j]]])];, {i,
Np}, {j, Np}];
Table[C1[[i,
i]] = -Total[Cases[C1[[i, All]], Except[C1[[i, i]]]]];, {i, Np}];
W1 = C1; W1[[1, All]] = 0.;
W2 = C1 . W1; W2[[Np, All]] = 0.;
W3 = C1 . W2; W3[[Np, All]] = 0.;
W4 = C1 . W3; W4[[1, All]] = 0.;
C0 = 1.8751^-4; K1M = C0 W4;
Y1[t_] := Table[Subscript[x, i][t], {i, 1, Np}]; Y4 = K1M . Y1[t];
\[Alpha] = -0.001; Xb =
2 G1 Cos[2 \[CapitalOmega] t] (Subscript[x, Np][
t] + \[Alpha] Subscript[x, Np][t]^3); force = (D[Xb, t,
t] + \[Mu] D[Xb, t]) Table[KroneckerDelta[1, i], {i, 1, Np}];
eqns = Thread[D[Y1[t], t, t] + \[Mu] D[Y1[t], t] + Y4 == force]; eq1 =
eqns[[1]] /.
Solve[Last[eqns], (Subscript[x, 10]^\[Prime]\[Prime])[t]]; eqns1 =
Join[{eq1}, Rest[eqns]];
Lg = 1; bt = 1.8751/Lg; \[Xi] = X;
y0 = -0.5 a (((Cos[bt*\[Xi]] - Cosh[bt*\[Xi]]) -
0.734*(Sin[bt*\[Xi]] - Sinh[bt*\[Xi]]))); X0 = -y0; X0d = 0 y0;
s = NDSolve[{eqns1, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax}];
Table[Plot[Evaluate[Subscript[x, i][t] /. First@s], {t, 0, tmax},
PlotRange -> All, PlotLabel -> Row[{"i = ", i}]], {i, 1, Np}]
Xbऔर$x_1(t)$लगभग 0.3 के अंतर से भिन्न होते हैं। हम इस स्थिरांक को प्रारंभिक स्थिति में शामिल कर सकते हैं$x_1(t)$ लेकिन साथ रहना बेहतर हो सकता है $x_1(0)=0$जैसा कि ऊपर दिए गए कोड में है। हमें यह भी ध्यान देना चाहिए, कि प्रस्तावित एल्गोरिथ्म मनमानी के लिए स्थिर नहीं है Np, लेकिन हम इसे बढ़ाकर स्थिर बना सकते हैं$\mu$ सीमा बिंदु के लिए $x_1$ जैसा कि आमतौर पर हमने लाइनों की विधि में किया था।
{Plot[{Evaluate[Xb /. First@s /. t -> t],
Evaluate[Subscript[x, 1][t] /. First@s]}, {t, 0, tmax},
PlotRange -> All],
Plot[{Evaluate[Xb /. First@s /. t -> t] -
Evaluate[Subscript[x, 1][t] /. First@s]}, {t, 0, tmax},
PlotRange -> All]}
thisstep = 0;
laststep = 0;
stepsize = 0;
First@NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax},
MaxStepFraction -> 1/15,
StepMonitor :> (laststep = thisstep; thisstep = t;
stepsize = thisstep - laststep;),
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic},
Method -> {"EventLocator",
"Event" :> (If[stepsize < 10^-9, 0, 1])}},
WorkingPrecision -> 24.05]
ReImPlot[#, {t, 0, laststep},
PlotRange -> {{0, laststep}, {900, -900}},
ColorFunction -> "DarkRainbow",
PlotLabels -> {"Expressions", "ReIm"}] & /@ %183[[All, 2]]
laststep
(* 7.12986 *)
ReImPlot[#, {t, 0, 7}, PlotRange -> Full,
ColorFunction -> "DarkRainbow"] & /@ %183[[2 ;; 9, 2]]
ReImPlot[%183[[1, 2]], {t, 0, laststep}, PlotRange -> Full,
ColorFunction -> "DarkRainbow",
PlotLabels -> {"Expressions", "ReIm"}]
StepDataPlot[%183]
पहला चैनल सबसे तेजी से दोलन करता है और आयाम तेजी से बढ़ रहा है। लक्ष्यों या परिशुद्धता के लिए प्रत्येक विधि विकल्प भारी शक्ति के साथ दोलनों की गणना करने में सक्षम है ताकि अन्य सभी चैनल केवल तेजी से बढ़ें। एक श्रेणी के रूप में इस स्थिरांक की गणना में दोलन होते हैं।
अनुकूलन समाधान के लिए सबसे लंबे डोमेन के परिप्रेक्ष्य से किया जाता है। चूंकि सभी समाधान चैनलों का वर्चस्व है$x_{1}$ यह सबसे महत्वपूर्ण है।
डोमेन को काटना जानकारीपूर्ण दृश्य के लिए अनुमति देता है:
ReImPlot[%183[[1, 2]], {t, 0, 4.3}, PlotRange -> Full,
ColorFunction -> "DarkRainbow",
PlotLabels -> {"Expressions", "ReIm"}]
का हल $x_{1}$समय-निर्भर आवृत्ति के साथ धीमी दोलन से मिलकर बनता है कि आवृत्ति समय के साथ तेज हो जाती है। यह धीरे-धीरे समय के साथ धीमी गति से लेकिन बहुत तेज आवृत्ति के साथ इस धीमे आवरण के नीचे दोलन करता है।
कथानक में शोर के लिए अनुमान के कारण कथानक को प्रभावित किया जाता है। ColorFunctionशो दोलन शून्य के माध्यम से जाना है। लिफाफा एक्स-अक्ष के संबंध में आयामों में असममित है।
यह एक मौका है कि 7.12986 में एकवचन और थोड़ी देर बाद एक बढ़ाया कार्यप्रणाली के साथ स्थिर गणना की जा सकती है।
सबसे अच्छा दृष्टिकोण हैं
sol = NDSolve[{eqns, Y1[0] == X0, Y1'[0] == X0d}, Y1[t], {t, 0, tmax},
Method -> {"Extrapolation", Method -> "ExplicitModifiedMidpoint",
"StiffnessTest" -> False}, WorkingPrecision -> 32];
ReImPlot[%198[[1, 1, 2]], {t, 0, 4.3}, PlotRange -> Full,
ColorFunction -> "DarkRainbow"]
ReImPlot[#, {t, 0, 7}, PlotRange -> Full,
ColorFunction -> "DarkRainbow"] & /@ %198[[1, 2 ;; 9, 2]]
दोनों विधियों के बीच, केवल थोड़ा अंतर है दोनों उच्च-परिशुद्धता हैं। लेकिन उपाय अलग हैं। दोनों तेजी से दोलनों से शोर और त्रुटि की गणना करते हैं। लेकिन समय में समाधान छोटे कदम और अधिक त्रुटि और शोर को जोड़ते हैं।
एक्सट्रपलेशन महत्वपूर्ण समय पर बहुत तेजी से 7.12986डायवर्ज करता है लेकिन उप-केंद्रों में, अन्य चैनलों में समाधान कम भारी दोलन हैं। डोमेन के एक उपखंड में कम संचय के कारण कम संचय हो सकता है, माना जाता है कि सीमाओं को सावधानी के साथ लिया जाता है। एक्सट्रपलेशन को अपनाने से दोलन को सुचारू करके कम शोर और त्रुटि को एकीकृत करने का मौका है।
मेरी समस्या यह है कि NDSolve के लिए "Extrapolation" विधि उदाहरणों में अधूरी है। दूसरी तरफ गणितज्ञ आंतरिक रूप से बहुत कुछ करता है। वह भी दोनों प्रस्तुत तरीकों समूहों के बीच समानता की उच्च डिग्री के कारण हो सकता है। गणना बहुत तेज है। एक इष्टतम है WorkingPrecision। जिसे आगे नहीं बढ़ाया जा सकता है। डोमेन की लंबाई का एक इष्टतम मूल्य है। इससे मुझे संदेह होता है।
मुझे यह अवधारणा मिली है कि यह केवल परिमित ऊँचाई की एक नाड़ी है और यह वक्र शक्ति के विलोपन की प्रक्रिया को शांत करता है। आगे कोई अनर्थकारी घटना नहीं है। लेकिन विचलन बहुत छोटे चरणों में बहुत तेज है। कंप्यूटिंग हमेशा समाप्त होता है संदेश कठोरता के संदेश के समान होता है चरण आकार बहुत छोटा हो जाता है। अनुचित कठोरता स्विचन के परिहार के साथ इसे दूर नहीं किया जा सकता है।
सभी छोटे-समय के दोलनों के समुचित एकीकरण को इस उत्तर के लिए मेरे द्वारा प्रस्तुत किए जाने की तुलना में बहुत अधिक समय और कंप्यूटिंग शक्ति की आवश्यकता हो सकती है।
गणना किए गए डोमेन के पहले भाग में लाभ की अच्छी तरह से कल्पना की गई है:
अतिरिक्त रेखीय समाधान अधिक रेखीय उपप्रकार में बहुत कम दोलन कर रहे हैं। यह एक ही दोलनों को बहुत शुरुआत में और इससे अधिक उपनल में हासिल करता है$⁄pi$। Stiffnessswitching की तुलना में डोमेन की ऊपरी सीमा पर दोलन गति बहुत अधिक है। यह उस समाधान की तुलना है जिसे प्रश्न में चुना गया है।
StepDataPlotइन उपशाखाओं में कठोरता का पता चलता है। Inbetween कोई कठोरता का निष्पादन किया जाता है। यह इन संगणनाओं को प्रश्न की तुलना में बहुत तेज बनाता है।
रणनीतिक सवाल यह है कि क्या आयाम पर दोलनों $-30$त्रुटि / शोर माना जाता है या समाधान का हिस्सा हैं। चूंकि रन-कुट्टा विधि त्रुटि को शून्य करने के लिए डिज़ाइन की गई है और शोर बल्कि सवाल महत्वपूर्ण है। परिणाम इस विचार में बदलने योग्य है कि उप-अंतराल पर कंप्यूटिंग पूर्ण अंतराल पर कंप्यूटिंग के लिए एक अनुकूलन है।
NDSolve इस तरह के कैंसिलेशन को कुछ हद तक आंतरिक रूप से रद्द करता है। साहित्य में, इस घटना को विचलन के साथ अराजकता में इंद्रधनुष या पथ कहा जा सकता है। जैसा कि कदम के प्रोग्राम इवेंट कंट्रोल से लिया जा सकता है यह दृष्टिकोण काम नहीं करता है। यह एक सवाल से अनुकूलित है जहां NDSolve कई शाखाओं के साथ समाधान पर काम कर रहा है। इसने शाखाओं का पता नहीं लगाया।
उपखंड शायद सबसे अच्छा है, खासकर अगर आयाम बिल्कुल बड़ा है $15$। संख्यात्मक प्रयोगों के लिए मान प्रश्न से लिया गया है। सबसे अधिक संभावना ब्याज की अधिक है।
कुछ शोध आयोजित करने के लिए यह वास्तव में NDSolve के लिए विधि की समझ पर क्या कर रहा है :
[Flatten [Trace [NDSolve [system], TraceInternal -> True]] का चयन करें,! FreeQ [#, विधि | NDSolve`MethodData] और]
अपने आप से पूछें: वुल्फ्राम मैथेमेटिका एनडीएसोल्व फ़ंक्शन फ़ंक्शन क्या हैं?
"एडम्स" - पूर्वसूचक - 12 के माध्यम से आदेश 1 के साथ सुधारक एडम्स विधि "बीडीएफ" -
5 के माध्यम से आदेश 1 के साथ गियर पिछड़े विभेदन सूत्र। ImplicitRungeKutta "- मनमानी के परिवार - आदेश निहित रन - कुट्टा तरीके" SymplecticPartitionedRungeKutta "- interleaved Runge - अलग हैमिल्टनियन सिस्टम के लिए कुट्टा विधियाँ" MethodOfLines "- PDEs" एक्सट्रापोलिशन "-" ग्रैप "के समाधान के लिए लाइनों की विधि। , संभव submethods के साथ [बुलेट] "ExplicitEuler" - फॉरवर्ड यूलर विधि [Bullet] "ExplicitMidpoint" - मिडपॉइंट नियम विधि [Bullet] "ExplicitModifiedMidpoint" - ग्रैग स्मूथिंग के लिए मिडपॉइंट नियम विधि [Bullet] "LinearlyImplicitEuler" - रैखिक रूप से निहित यूलर विधि ) "डबलस्टेप" - "एक्सट्रैपलेशन" "लोकेएक्सएक्ट" का "बेबी" संस्करण - स्थानीय रूप से सटीक प्रतीकात्मक समाधान "स्टिफनेसस्विचिंग" के लिए संख्यात्मक सन्निकटन - एकीकरण के बीच में नॉनस्टिफ और कठोर तरीकों के बीच स्विच करने की अनुमति देता है
"प्रोजेक्शन" - अपरिवर्तनीय - संरक्षण विधि "ऑर्थोगोनलप्रोजेक्शन "- विधि जो समाधानों की अलंकारिकता को संरक्षित करती है" आईडीए "-
अंतर के सिस्टम के लिए प्रारंभिक मूल्य की समस्या के लिए सामान्य उद्देश्य सॉल्वर - बीजीय समीकरण (डीएई)" शूटिंग "- बीवीपी के लिए शूटिंग विधि" पीछा करना "- गेलैंड - बीवीपी के लिए लोकसिटीवस्कैसी पीछा करने की विधि" EventLocator "- असंतोष, अवधियों, आदि का पता लगाने के लिए घटना का स्थान" परिमितियां "- परिमित तत्व समस्याएं
निगरानी और एल्गोरिदम का चयन करें :
try[m_] :=
Block[{s = e = 0},
NDSolve[system, StepMonitor :> s++, EvaluationMonitor :> e++,
Method -> m]; {s, e}]
वास्तविक ब्याज और अच्छे समाधान की विधि और विकल्प के साथ। यह दुखद है कि यह ट्यूटोरियल वास्तव में गहराई तक जाता है। चयन प्रक्रिया में काफी समय लगता है।
यह प्रदर्शन पक्ष की कार्यप्रणाली को दर्शाता है: एक 3D सतह की साजिश रचने के कार्य में अनुकूली 3D प्लॉटिंग । यह स्वयं स्टीफन वोल्फ्राम का एक प्रदर्शन है। और इस के अधिक हैं। एक्स-प्लॉटिंग के लिए एक है: अनुकूली प्लॉटिंग । यह ट्यूटोरियल NDSolve के लिए "DoubleStep" विधि दिखाता है । यह इस समस्या के लिए रन-कुट्टा-विधि सफल क्यों है, इस पर एक नज़र देता है। यह ट्यूटोरियल कुछ हद तक किसी भी तरह से पाठक को उस Methodविकल्प के पीछे छिपे हुए कॉम्प्लेक्स "Automatic"में ले जाता है, जो साहित्य, गणितज्ञ प्रलेखन में NDSolve समाधान में इतना सर्वव्यापी है। अच्छा अभ्यास यह है कि प्लॉट फ़ंक्शन के रूप में अनुकूली नमूना कैसे प्राप्त किया जाए ।
समस्या यह है कि "NIntegrate के लिए आपको सांख्यिक मूल्यांकन को बाध्य करना चाहिए" के द्वारा चिह्नित के समान है, अन्यथा यह कुछ क्वाड्रेचर स्कीम को नियोजित कर सकता है जो eval अंक की संख्या को कम करता है। "
तथा
"एकीकृत NIntegrate के प्रतीकात्मक रूप से मूल्यांकन बिंदुओं की संख्या को कम करने के लिए आवधिकता जैसी कुछ विशेषताओं का पता लगा सकता है। AFAIK यह प्रतीकात्मक लागू करेगा जब तक कि आप इसे विधि से बंद नहीं करते हैं -> {स्वचालित," प्रतीक चिह्न " (स्वत: के बजाय स्पष्ट विधि विनिर्देश हो सकता है) या "ब्लैक बॉक्स" विधि (_? NumericQ) का उपयोग करके। दोनों तरीके चतुर्भुज योजना को अक्षम नहीं करते हैं।
उपखंड के लिए एक अच्छी अवधारणा इस समुदाय से 2d में कार्यों की गणना करने के लिए धीमी गति से समोच्च भूखंड अनुकूली नमूने को गति देने के लिए दी गई है ।
दिए गए डेटा के साथ दी गई समस्या ग़ैर-कठोर है, लेकिन कठोर हो जाती है यदि प्रश्न से सटीक विकल्प उस कठोर को लिया जाता है। जैसा कि गणितज्ञ प्रलेखन का अध्ययन करके पुष्टि की जा सकती है कि सिफारिश का विकल्प पूरी तरह से है WorkingPrecision।
इंटरपोलटिंगफंक्शन के कई उदाहरणों को एक साथ विभाजित करने के तरीके के साथ काम करें ! आगे का महत्वपूर्ण चरण प्रत्येक पूर्ण अवधि को ठीक से ध्यान में रखना है। उपखंड के लिए अच्छा तरीका NDSolve EventLocator में प्रलेखित है