Metode kuadratur diferensial gagal pada PDE orde 4 dengan bc nonlinier karena grid semakin padat

Dec 26 2020

Saya menggunakan metode kuadratur diferensial (DQM) untuk memecahkan masalah nilai batas awal berikut:

$$\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$$

Sini $X_b(t)$( Xb[t]dalam kode di bawah) adalah input dalam sistem. Untuk fungsi harmonik

$$X_b(t)=2G_1 \cos(2\Omega t)$$

sebagai masukan, NDSolvebekerja dengan sempurna. Untuk input lain juga, simulasi berjalan dengan baik. Tapi untuk masukannya

$$X_b(t)=2G \cos(2 \Omega t) (w(1,t)+\alpha \; w(1,t)^3)$$

Osilasi frekuensi tinggi tumbuh dan hasil simulasi menjadi semakin tidak akurat, karena jumlah titik kisi ( Npdalam kode di bawah) meningkat. Jika tmaxbesar atau Np > 10, maka simulasi gagal dengan peringatan

NDSolve :: ndsz: diduga sistem kaku atau singularitas.

Saya telah memecahkan masalah ini dengan metode lain, tidak ada osilasi frekuensi tinggi.

Berikut ini adalah percobaan saya. PDE telah Np - 1dialihkan ke ODE dengan DQM.

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;

Input Xb[t]diganti dalam sistem persamaan melalui vektor kolom 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]

Hasil yang diperoleh dengan versi 11.1 untuk Np = 6dan 8diberikan di bawah ini. Sebab Np = 8, fluktuasi output meningkat.

Sebab Np = 12, itu memberi peringatan

NDSolve :: ndsz: Pada t == 7.129860412016928`, ukuran langkah secara efektif nol; singularitas atau sistem kaku yang dicurigai.

Saya telah menggunakan opsi yang berbeda NDSolveuntuk menangani sistem yang kaku, tetapi tetap tidak berfungsi.

Saya pikir, jika saya gunakan NDSolvedalam interval yang lebih kecil maka mungkin berhasil. Jadi saya membuat kode di mana kondisi awal (untuk iterasi ke-i) dihitung berdasarkan keluaran dari iterasi sebelumnya (i - 1). Tetapi beberapa NDSolvesimulasi juga tidak berhasil.

Saya telah mencoba kondisi awal yang berbeda juga, tetapi peringatannya tetap ada. Adakah yang bisa membantu saya untuk memecahkan masalah ini? Terima kasih.

Jawaban

9 xzczd Jan 03 2021 at 11:40

Jika DQM diimplementasikan dengan benar, maka ini mungkin menjadi batasan penting dari metode ini. Saya tidak tahu apa-apa tentang DQM, tetapi memindai makalah ini , saya merasa metodenya mirip Pseudospectral. Memang, uji cepat menunjukkan bahwa, matriks koefisien pembobotan turunan orde 1 di DQM sama persis dengan matriks diferensiasi turunan orde 1 dalam metode pseudospektral:

NDSolve`FiniteDifferenceDerivative[1, X, DifferenceOrder -> "Pseudospectral"][
  "DifferentiationMatrix"] == C1[[All, All, 1]]
(* True *)

Meskipun saya tidak dapat memberikan contoh konkret saat ini, saya mengamati bahwa hal itu Pseudospectralmungkin menjadi tidak stabil ketika titik kisi spasial meningkat dalam kasus tertentu. Mari kita uji apakah masalah Anda termasuk dalam hal semacam itu. Karena bc khusus, kita tidak dapat menggunakan NDSolvesecara langsung untuk menyelesaikan masalah, jadi mari diskritkan sistem di$x$arah sendiri. Saya akan gunakan pdetoodeuntuk tugas itu.

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: Pada t == 2.356194489774355`, ukuran langkah efektif nol; singularitas atau sistem kaku yang dicurigai.

Seperti yang bisa kita lihat, ketika jumlah titik grid tidak lebih dari 10, solusinya tampaknya stabil dan konvergen cukup cepat, tetapi begitu pointsmeningkat menjadi 11, solusi menjadi liar, mirip dengan perilaku implementasi OP.

Lantas, bagaimana cara mengelak? Menggunakan rumus perbedaan orde rendah untuk mendiskrit ternyata efektif:

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

Seperti yang ditunjukkan di atas, solusinya tetap stabil dengan points = 50; difforder = 2.

difforder = 4 juga dapat digunakan jika Anda suka.


Tambahan: Penerapan kembali Metode OP

Setelah melihat lebih dekat kode OP bersama dengan kertas yang ditautkan di awal jawaban ini, saya rasa saya sudah mengerti apa yang telah diterapkan OP. Berikut ini adalah implementasi saya untuk metode yang sama, yang mungkin sedikit lebih mudah dibaca:

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]

Beberapa penjelasan lebih lanjut: dalam metode ini file $\frac{\partial^4}{\partial x^4}$ telah didiskritisasi dengan cara rekursif yaitu $\frac{\partial}{\partial x}$didiskritisasi terlebih dahulu ( C1[[All, All, 1]]dalam kode OP, d1dalam kode saya) dan diskritisasi$\frac{\partial^4}{\partial x^4}$dihitung menggunakan Dot. Masih merasa curiga? Oke, mari kita validasi:

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'''"}]

Sejak $\frac{\partial}{\partial x}$, $\frac{\partial^2}{\partial x^2}$ dan $\frac{\partial^3}{\partial x^3}$ semua muncul sebagai perantara dalam metode ini, masalah bcs OP dapat diterapkan dengan memodifikasi matriks secara langsung, misalnya:

ListLinePlot[{grid, d1withbc . f[grid]}\[Transpose]]

Seperti yang diilustrasikan di atas, $\frac{\partial f}{\partial x}\Big|_{x=0}=0$ telah diberlakukan.

7 AlexTrounev Dec 28 2020 at 06:04

Karena kode ini merupakan implementasi DQM untuk balok kantilever maka kita perlu menempatkan kondisi batas yang tepat agar kode ini stabil dengan jumlah titik grid yang Npberubah. Ini hanya modifikasi kecil tetapi berfungsi untuk semua Np, misalnya

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

Dengan pendekatan ini kita harus mempertimbangkan Xb[t]sebagai kekuatan eksternal yang diterapkan pada titik arbitrer nextsebagai

Xb[t] = 2 G1 Cos[
       2 Ω t] (Subscript[x, Np][
         t] + α Subscript[x, Np][t]^3) Table[
       KroneckerDelta[next, i], {i, 2, Np}]; 

Dalam kasus next=Npkami memiliki kode di atas. Alasan utama mengapa kode asli menghasilkan solusi yang tidak stabil adalah K1Mdefinisi matriks yang diambil dari makalah Application of the Generalized Differential Quadrature Method to the Study of Pull-In Phenomena dari MEMS Switches, oleh Hamed Sadeghian, Ghader Rezazadeh, dan Peter M. Osterberg. Kita dapat mendefinisikannya kembali sebagai berikut

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}];

Kita dapat membandingkan solusi ini di Xb=0(titik merah) dengan solusi yang dihasilkan oleh kode xzczd dengan points=10(garis penuh)

Nah jika kita meletakkan Np=30dan mengaplikasikannya Xbke titik pertama seperti pada kode di atas, maka kita mendapatkan gambar untuk setiap titik grid sebagai berikut

Table[Plot[Evaluate[Subscript[x, i][t] /. First@s], {t, 0, tmax}, 
  PlotRange -> All, PlotLabel -> Row[{"i = ", i}]], {i, 1, Np}]

Ini sangat umum menanggapi gaya eksternal. Dengan menggunakan matriks ini K1M = C0 W4kita dapat mewujudkan gagasan utama Xbimplementasi sebagai$x_1(t)$ sebagai berikut

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

Akhirnya kita bisa memeriksanya Xbdan$x_1(t)$berbeda dengan konstanta sekitar 0,3. Kita bisa memasukkan konstanta ini pada kondisi awal untuk$x_1(t)$ tapi bisa lebih baik untuk tetap bersama $x_1(0)=0$seperti pada kode diatas. Juga harus kita catat, bahwa algoritma yang diusulkan tidak stabil untuk sembarang Np, tapi kita bisa membuatnya stabil dengan meningkatkannya$\mu$ untuk titik batas $x_1$ seperti yang biasa kami lakukan dalam metode garis.

{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]}  

2 SteffenJaeschke Jan 04 2021 at 04:41
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]

Saluran pertama berosilasi paling cepat dan amplitudo membesar secara eksponensial. Setiap opsi metode untuk tujuan atau presisi mampu menghitung osilasi dengan kekuatan yang sangat besar sehingga semua saluran lainnya hanya tumbuh secara eksponensial. Dalam bentuk rentang yang menghitung konstanta ini ada osilasi.

Optimasi dilakukan dengan perspektif domain terpanjang untuk solusinya. Karena semua saluran solusi didominasi oleh$x_{1}$ itu yang paling penting.

Memotong domain memungkinkan tampilan informatif:

ReImPlot[%183[[1, 2]], {t, 0, 4.3}, PlotRange -> Full, 
 ColorFunction -> "DarkRainbow", 
 PlotLabels -> {"Expressions", "ReIm"}]

Solusi dari $x_{1}$terdiri dari osilasi yang lebih lambat dengan frekuensi yang bergantung pada waktu sehingga frekuensi semakin cepat seiring waktu. Ini berosilasi di bawah selubung yang lebih lambat ini dengan lambat dengan waktu yang semakin berkurang tetapi frekuensi yang jauh lebih cepat.

Plot tidak tepat karena adanya kesimpulan hingga kebisingan dalam plot. Ini ColorFunctionmenunjukkan bahwa osilasi melewati nol. Amplopnya asimetris dalam amplitudo terhadap sumbu x.

Ini adalah kemungkinan bahwa singularitas di 7.12986 dan beberapa saat kemudian dapat dihitung stabil dengan metodologi yang disempurnakan.

Pendekatan terbaik adalah

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]]

Di antara kedua metode tersebut, hanya ada sedikit perbedaan keduanya adalah presisi tinggi. Tetapi solusinya berbeda. Baik menghitung noise dan error paling banyak dari osilasi cepat. Tetapi semakin kecil solusi langkah dalam waktu semakin banyak kesalahan dan kebisingan bertambah.

Ekstrapolasi menyimpang jauh lebih cepat pada waktu kritis 7.12986tetapi dalam sub-interval, solusi di saluran lain tidak terlalu berosilasi. Sebuah subdivisi dari domain dapat menyebabkan osilasi yang lebih sedikit karena lebih sedikit akumulasi lentur yang seharusnya dilakukan dengan hati-hati. Ada peluang untuk mengintegrasikan lebih sedikit noise dan error dengan menghaluskan osilasi dengan mengadopsi ekstrapolasi.

Masalah saya adalah bahwa Metode "Ekstrapolasi" untuk NDSolve tidak lengkap dalam contoh. Mathematica di sisi lain melakukan banyak hal secara internal. Itu mungkin juga karena tingkat kemiripan yang tinggi antara kedua kelompok metode yang disajikan. Perhitungannya sangat cepat. Ada yang optimal WorkingPrecision. Itu tidak bisa ditingkatkan lebih jauh. Panjang domain memiliki nilai yang optimal. Itu membuat saya skeptis.

Saya telah mendapatkan konsep yang hanyalah denyut nadi dengan ketinggian yang terbatas dan bahwa kurva tersebut menjadi tenang dalam proses pemusnahan listrik. Tidak ada peristiwa bencana di depan. Tetapi perbedaannya sangat cepat, banyak urutan dalam langkah-langkah yang sangat kecil. Komputasi selalu berakhir adalah pesan yang mirip dengan pesan kekakuan ukuran langkah menjadi terlalu kecil. Itu tidak dapat diatasi dengan menghindari peralihan kekakuan yang tidak tepat.

Integrasi yang tepat dari semua osilasi kecil-kecilan mungkin memerlukan lebih banyak waktu dan daya komputasi daripada yang dapat saya tawarkan untuk jawaban ini.

Keuntungan di bagian pertama dari domain yang dihitung divisualisasikan dengan baik oleh:

Solusi yang diekstrapolasi jauh lebih sedikit berosilasi dalam subinterval yang lebih linier. Ini memperoleh osilasi yang sama di awal dan di subinterval lebih besar dari$⁄pi$. Momentum osilasi jauh lebih tinggi di batas atas domain dibandingkan dengan Stiffnessswitching. Ini adalah perbandingan solusi yang dipilih dalam pertanyaan.

Mengevaluasi StepDataPlotmenunjukkan bahwa dalam sub-interval ini pengalihan kekakuan terjadi. Antara tidak ada kekakuan yang dieksekusi. Ini membuat penghitungan ini jauh lebih cepat daripada penghitungan dari pertanyaan.

Pertanyaan strategisnya adalah apakah osilasi pada amplitudo $-30$dianggap error / noise atau merupakan bagian dari solusi. Karena metode Runge-Kutta dirancang untuk memusatkan perhatian pada kesalahan dan kebisingan, pertanyaannya agak penting. Hasilnya dapat diubah menjadi gagasan bahwa komputasi pada subinterval adalah pengoptimalan untuk menghitung selama interval lengkap.

NDSolve melakukan pembatalan tersebut secara internal sampai batas tertentu. Dalam literatur, fenomena ini bisa disebut pelangi atau jalan menuju kekacauan dengan divergensi. Seperti yang dapat diambil dari kontrol acara yang diprogram dari ukuran anak tangga, pendekatan ini tidak berfungsi. Ini diadaptasi dari pertanyaan di mana NDSolve beroperasi pada solusi dengan banyak cabang. Itu tidak mendeteksi cabang sama sekali.

Pembagian ini mungkin yang terbaik terutama jika amplitudo benar-benar lebih besar dari $15$. Nilai untuk eksperimen numerik diambil dari pertanyaan. Kemungkinan besar ada lebih banyak minat.

Untuk melakukan beberapa penelitian apa yang sebenarnya dilakukan, lihat pemahaman metode untuk NDSolve :

Pilih [Ratakan [Trace [NDSolve [system], TraceInternal -> True]],! FreeQ [#, Metode | NDSolve`MethodData] &]

Tanyakan pada Diri Anda: Apa saja metode fungsi Wolfram Mathematica NDSolve?

"Adams" - prediktor - korektor Metode Adams dengan urutan 1 hingga 12 "BDF" - Persiapkan rumus diferensiasi mundur implisit dengan urutan 1
hingga 5 "ExplicitRungeKutta" - pasangan tersemat adaptif dari 2 (1) hingga 9 (8) metode Runge - Kutta " ImplicitRungeKutta "- kelompok sembarang - urutan implisit Runge - metode Kutta" SymplecticPartitionedRungeKutta "- metode Runge - Kutta berseling untuk sistem Hamiltonian yang dapat dipisahkan" MethodOfLines "- metode garis untuk solusi PDE" Ekstrapolasi "- (Gragg -) Bulirsch - Metode ekstrapolasi Stoer , dengan kemungkinan submetode [Bullet] "ExplicitEuler" - metode Euler maju [Bullet] "ExplicitMidpoint" - metode aturan titik tengah [Bullet] "ExplicitModifiedMidpoint" - metode aturan titik tengah dengan penghalusan Gragg [Bullet] "LinearlyImplicitEuler" - metode Euler implisit linier [Bullet ] "LinearlyImplicitMidpoint" - metode aturan titik tengah implisit linier [Bullet] "LinearlyImplicitModifiedMidpoint" - Bader implisit linier - metode aturan titik tengah yang dihaluskan "DoubleStep" - versi "baby" dari "Extrapolation" "LocallyExact" - perkiraan numerik untuk solusi simbolik yang tepat secara lokal "StiffnessSwitching" - memungkinkan peralihan antara metode nonstiff dan stiff di tengah
integrasi "Projection" - invariant - preserving method "OrthogonalProjection "- metode yang mempertahankan orthonormalitas solusi" IDA "- pemecah tujuan umum untuk masalah nilai awal untuk
sistem persamaan aljabar - diferensial (DAE)" Menembak "- metode pemotretan untuk BVP" Mengejar "- Gelfand - Metode pengejaran Lokutsiyevskii untuk BVP" EventLocator "- lokasi acara untuk mendeteksi diskontinuitas, titik, dll." FiniteElements "- masalah elemen hingga

Gunakan Pemantauan dan Pemilihan Algoritma :

try[m_] := 
 Block[{s = e = 0}, 
  NDSolve[system, StepMonitor :> s++, EvaluationMonitor :> e++, 
   Method -> m]; {s, e}]

dengan Metode dan pilihan minat nyata dan solusi yang baik. Sayangnya tutorial ini benar-benar mendalam. Proses seleksi memakan banyak waktu.

Demonstrasi ini menunjukkan metodologi yang disukai: Plotting 3D Adaptif pada tugas merencanakan permukaan 3D. Ini adalah demonstrasi dari Stephen Wolfram sendiri. Dan masih ada lagi ini. Ada satu untuk xy-plotting: Adaptive Plotting . Tutorial ini menunjukkan Metode "DoubleStep" untuk NDSolve . Ini menawarkan pandangan mengapa metode Runge-Kutta berhasil untuk masalah ini. Tutorial ini agak dan entah bagaimana mendorong pembaca ke kompleks yang tersembunyi di balik Methodopsi "Automatic"yang ada di mana-mana dalam solusi NDSolve dalam literatur, dokumentasi Mathematica. Praktik yang baik adalah bagaimana mendapatkan sampling adaptif seperti pada fungsi plot .

Masalahnya mirip dengan yang dilambangkan dengan "untuk NIntegrate Anda harus memaksa evaluasi numerik, jika tidak, mungkin menggunakan beberapa skema kuadrat yang meminimalkan jumlah poin eval."

dan

"Dari bentuk simbolik integrand NIntegrate dapat mendeteksi beberapa fiturnya seperti periodisitas untuk meminimalkan jumlah titik evaluasi. AFAIK akan menerapkan simbolik sampai Anda mematikannya dengan Metode -> {Otomatis," SymbolicProcessing "-> Tidak Ada} (Alih-alih Otomatis mungkin merupakan spesifikasi metode eksplisit) atau dengan menggunakan metode "kotak hitam" (_? NumericQ). Kedua cara tersebut tidak menonaktifkan skema kuadrat. "

Sebuah konsep bagus untuk subdivisi diberikan dalam percepatan plot kontur pengambilan sampel adaptif untuk fungsi komputasi lambat dalam 2d dari komunitas ini.

Masalah yang diberikan dengan data yang diberikan tidak kaku tetapi menjadi kaku jika opsi presisi dari pertanyaan diambil yang kaku. Seperti yang dapat dikonfirmasi dengan mempelajari dokumentasi Mathematica, pilihan rekomendasi semata-mata WorkingPrecision.

Bekerja dengan cara menggabungkan beberapa contoh fungsi interpolasi ! Langkah penting ke depan adalah setiap periode penuh harus diperhitungkan dengan benar. Metode bagus untuk subdivisi didokumentasikan di NDSolve EventLocator