グリッドが密になるにつれて、微分求積法は非線形bcの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：特異性またはスティッフなシステムが疑われます。
私は別の方法でこの問題を解決しました。高周波振動は存在しません。
以下は私の試用版です。PDEはNp - 1、DQMを使用してODEに離散化されています。
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 = 68Np = 8
の場合Np = 12、警告が表示されます
NDSolve :: ndsz：t == 7.129860412016928`では、ステップサイズは事実上ゼロです。特異性または硬いシステムが疑われる。
NDSolve硬いシステムに対処するためにさまざまなオプションを使用しましたが、それでも機能しません。
NDSolve間隔を短くすれば成功するのではないかと思いました。そこで、前の反復（i-1）からの出力に基づいて初期条件（i番目の反復）を計算するコードを作成しました。しかし、複数のNDSolveシミュレーションも機能しませんでした。
さまざまな初期条件も試しましたが、警告が残ります。誰かが私が問題を解決するのを手伝ってくれますか？ありがとう。
回答
DQMが正しく実装されている場合、これはメソッドの本質的な制限である可能性があります。DQMについては何も知りませんでしたが、この論文をスキャンすると、方法はに似ているように感じPseudospectralます。実際、簡単なテストでは、DQMの1次導関数の重み係数行列は、疑似スペクトル法の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：t == 2.356194489774355`では、ステップサイズは事実上ゼロです。特異性または硬いシステムが疑われる。
ご覧のとおり、グリッドポイントの数が10以下の場合、解は安定していて非常に速く収束しているように見えますが、にpoints増加する11と、OPの実装の動作と同様に、解はワイルドになります。
それで、どのように回避するのですか？低次の差の式を使用して離散化すると、効果的であることがわかります。
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 必要に応じて使用することもできます。
補遺：OPの方法の再実装
この回答の冒頭にリンクされているペーパーと一緒にOPのコードを詳しく調べた後、OPが何を実装したかを理解したと思います。以下は、同じメソッドの実装ですが、もう少し読みやすいかもしれません。
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]]OPのコード、d1私のコードで）、離散化されます$\frac{\partial^4}{\partial x^4}$を使用して計算されDotます。それでも疑わしいと感じますか？OK、検証しましょう：
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}$ すべてがメソッドの中間として表示されている場合、OPの問題のbcsは、マトリックスを直接変更することによって課すことができます。次に例を示します。
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、我々上記のコードを持っています。元のコードが不安定な解を生成する主な理由は、K1MHamed Sadeghian、Ghader Rezazadeh、およびPeter M.Osterbergによる論文「MEMSスイッチのプルイン現象の研究への一般化微分直交法の適用」から取られた行列定義です。次のように再定義できます
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（赤い点）でxzczdコードによって生成されたソリューションとpoints=10（実線）で比較できます。
上記のコードのように最初のポイントを配置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振動がゼロを通過することを示しています。エンベロープは、x軸に関して振幅が非対称です。
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の「外挿」方法が例では不完全であるということです。反対側のMathematicaは内部的に非常に多くのことをします。これも、提示された両方のメソッドグループ間の高度な類似性のためである可能性があります。計算は非常に高速です。最適がありWorkingPrecisionます。それ以上強化することはできません。ドメインの長さには最適な値があります。それは私を懐疑的にさせます。
私は、有限の高さのパルスであり、電力が消滅する過程で曲線が落ち着くという概念を持っています。壊滅的な出来事はありません。しかし、発散は非常に速く、非常に小さなステップで多くの順序になります。計算は常に終了しますが、ステップサイズが小さくなりすぎる剛性メッセージに似たメッセージです。これは、不適切な剛性切り替えを回避することでは克服できません。
すべての短時間の振動を適切に統合するには、この答えに提供できるよりもはるかに多くの時間と計算能力が必要になる場合があります。
計算された領域の最初の部分の利点は、次のようによく視覚化されます。
外挿された解は、より線形なサブインターバルでは振動がはるかに少なくなります。それは、最初とサブインターバルで同じ振動を獲得します。$⁄pi$。振動運動量は、剛性スイッチングよりもドメインの上限ではるかに高くなります。これは、質問で選択されたソリューションの比較です。
を評価するStepDataPlotと、これらのサブインターバルで剛性切り替えが行われることがわかります。その間、スティフネスウィッチングは実行されません。これにより、これらの計算は質問の計算よりもはるかに高速になります。
戦略的な問題は、振幅での振動かどうかです $-30$エラー/ノイズと見なされるか、ソリューションの一部です。ルンゲクッタ法はエラーとノイズをゼロにするように設計されているため、問題はかなり重要です。結果は、サブインターバルでの計算が完全なインターバルでのコンピューティングの最適化であるという考えに変換できます。
NDSolveは、すでにある程度内部でそのようなキャンセルを行っています。文献では、この現象は虹または発散を伴う混沌への道と呼ばれることがあります。ステップサイズのプログラムされたイベント制御から取得できるように、このアプローチは機能しません。これは、NDSolveが多くのブランチを持つソリューションで動作しているという質問から適応されています。ブランチはまったく検出されませんでした。
特に振幅が絶対に大きい場合は、細分割がおそらく最適です。 $15$。数値実験の値は質問から取得されます。最も可能性が高いのは、もっと興味深いものです。
これが実際に何をしているのかを調査するために、NDSolveの方法の理解を見てください。
Select [Flatten [Trace [NDSolve [system]、TraceInternal-> True]]、！FreeQ [＃、メソッド| NDSolve`MethodData]＆]
自問してみてください：Wolfram Mathematica NDSolve関数メソッドとは何ですか？
「アダムス」-予測子-修正子アダムス法（1次から12次）「BDF」-1
から5次のギア暗黙的後方微分式「ExplicitRungeKutta」-2（1）から9（8）ルンゲ-クッタ法の適応埋め込みペア " ImplicitRungeKutta "-任意のファミリー-順序暗黙的ルンゲ-Kuttaメソッド" SymplecticPartitionedRungeKutta "-インターリーブされたルンゲ-分離可能なハミルトニアンシステムのKuttaメソッド" MethodOfLines "-PDEの解法のラインメソッド" Extrapolation "-（Gragg-）Bulirsch-Stoer外挿法、可能なサブメソッドを使用[Bullet] "ExplicitEuler"-フォワードオイラーメソッド[Bullet] "ExplicitMidpoint"-中点ルールメソッド[Bullet] "ExplicitModifiedMidpoint"-グラッグスムージングを使用した中点ルールメソッド[Bullet] "LinearlyImplicitEuler"-線形暗黙オイラーメソッド[Bullet ] "LinearlyImplicitMidpoint"-線形に暗黙的な中点法法[箇条書き] "LinearlyImplicitModifiedMidpoint"-線形に暗黙的なBader-平滑化された中点法 「DoubleStep」-「Extrapolation」の「baby」バージョン「LocallyExact」-局所的に正確なシンボリックソリューション「StiffnessSwitching」の数値近似-
統合の途中で非スティッフメソッドとスティッフメソッドを切り替えることができます「Projection」-不変-保存メソッド「OrthogonalProjection "-解の
正則性を維持する方法" IDA "-微分システムの初期値問題の汎用ソルバー-代数方程式（DAE）"射撃 "-BVPの射撃法"追跡 "-Gelfand-BVPのLokutsiyevskii追跡法" EventLocator」-不連続性、期間などを検出するためのイベントの場所「FiniteElements」-有限要素問題
監視および選択アルゴリズムを使用する：
try[m_] :=
Block[{s = e = 0},
NDSolve[system, StepMonitor :> s++, EvaluationMonitor :> e++,
Method -> m]; {s, e}]
真に関心のある方法とオプション、そして優れたソリューションを備えています。このチュートリアルが本当に深く掘り下げられているのは悲しいことです。選択プロセスには多くの時間がかかります。
このデモンストレーションは、好意的な方法論を示しています。3Dサーフェスをプロットするタスクでの適応3Dプロット。これはスティーブン・ウルフラム自身によるデモンストレーションです。そして、これはもっとあります。xyプロット用に1つあります：適応プロット。このチュートリアルでは、NDSolveの「DoubleStep」メソッドを示します。ルンゲクッタ法がこの問題に対して成功する理由を説明します。このチュートリアルは、読者を、文献、MathematicaドキュメントのNDSolveソリューションに非常に遍在するMethodオプションの背後に隠された複雑なものにいくらかそしてどういうわけか"Automatic"駆り立てます。グッドプラクティスはある関数plotのように適応サンプリングを取得する方法。
この問題は、「NIntegrateの場合、数値評価を強制する必要があります。そうしないと、評価ポイントの数を最小限に抑える求積スキームを採用する可能性があります」で示される問題と似ています。
そして
「被積分関数のシンボリック形式から、NIntegrateは、評価ポイントの数を最小限に抑えるために、周期性などのいくつかの機能を検出できます。方法-> {自動、「SymbolicProcessing」->なし}でオフにするまで、シンボリックを適用します。 （自動の代わりに明示的なメソッド指定である可能性があります）または「ブラックボックス」メソッド（_？NumericQ）を使用します。どちらの方法でも直交スキームが無効になることはありません。」
サブディビジョンの優れた概念は、このコミュニティから2Dで関数を計算するのが遅い場合の等高線図の適応サンプリングの高速化に示されています。
与えられたデータに関する与えられた問題は堅くありませんが、質問からの精度オプションがその堅いものであると、固くなります。Mathematicaのドキュメントを調べることで確認できるように、推奨の選択はただWorkingPrecisionです。
仕事を一緒れるInterpolatingFunctionの複数のインスタンスをスプライスする方法！今後の重要なステップは、各全期間を適切に考慮に入れる必要があることです。細分化の優れた方法は、NDSolveEventLocatorに記載されています。