O método de quadratura diferencial falha no PDE de 4ª ordem com não linear bc conforme a grade fica mais densa
Estou usando o método de quadratura diferencial (DQM) para resolver o seguinte problema de valor limite inicial:
$$\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$$
Aqui $X_b(t)$( Xb[t]
no código abaixo) é a entrada no sistema. Para função harmônica
$$X_b(t)=2G_1 \cos(2\Omega t)$$
como entrada, NDSolve
funciona perfeitamente. Para outras entradas também, a simulação continua corretamente. Mas para a entrada
$$X_b(t)=2G \cos(2 \Omega t) (w(1,t)+\alpha \; w(1,t)^3)$$
As oscilações de alta frequência aumentam e os resultados da simulação se tornam cada vez mais imprecisos, conforme o número de pontos de grade ( Np
no código abaixo) aumenta. Se tmax
for grande ou Np > 10
, a simulação falha com aviso
NDSolve :: ndsz: suspeita de singularidade ou sistema rígido.
Resolvi este problema com outro método, não há oscilações de alta frequência presentes.
O seguinte é meu julgamento. O PDE foi discretizado para Np - 1
ODEs com 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;
A entrada Xb[t]
é substituída no sistema de equações por meio de um vetor de coluna 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]
Os resultados obtidos com a versão 11.1 de Np = 6
e 8
são dadas abaixo. Pois Np = 8
, a flutuação na produção aumenta.

Pois Np = 12
, dá um aviso
NDSolve :: ndsz: Em t == 7.129860412016928`, o tamanho do passo é efetivamente zero; singularidade ou sistema rígido suspeito.

Usei diferentes opções de NDSolve
para lidar com o sistema rígido, mas ainda não está funcionando.
Pensei que, se usar NDSolve
em intervalos menores, pode dar certo. Então eu fiz o código no qual as condições iniciais (para a iª iteração) são calculadas com base na saída da iteração anterior (i - 1). Mas a NDSolve
simulação múltipla também não funcionou.
Também tentei diferentes condições iniciais, mas o aviso permanece. Alguém pode me ajudar a resolver o problema? Obrigado.
Respostas
Se o DQM for implementado corretamente, essa pode ser uma limitação essencial do método. Eu não sabia nada sobre DQM, mas digitalizando este papel , tenho a sensação de que o método é semelhante Pseudospectral
. Na verdade, um teste rápido mostra que, a matriz do coeficiente de ponderação da derivada de 1ª ordem em DQM é exatamente a mesma que a matriz de diferenciação da derivada de 1ª ordem no método pseudoespectral:
NDSolve`FiniteDifferenceDerivative[1, X, DifferenceOrder -> "Pseudospectral"][
"DifferentiationMatrix"] == C1[[All, All, 1]]
(* True *)
Embora eu não possa dar um exemplo concreto no momento, observei que Pseudospectral
pode se tornar instável quando os pontos da grade espacial aumentam em certos casos. Vamos testar se o seu problema pertence a esse tipo de coisa. Por causa do bc especial, não podemos usar NDSolve
diretamente para resolver o problema, então vamos discretizar o sistema em$x$direção por nós mesmos. Vou usar pdetoodepara a tarefa.
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: Em t == 2.356194489774355`, o tamanho do passo é efetivamente zero; singularidade ou sistema rígido suspeito.

Como podemos ver, quando o número de pontos da grade não passa de 10, a solução parece ser estável e convergir bastante rápido, mas uma vez que points
aumenta para 11
, a solução torna-se selvagem, semelhante ao comportamento da implementação do OP.
Então, como contornar? Usar a fórmula de diferença de ordem baixa para discretizar acaba sendo eficaz:
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}]

Conforme mostrado acima, a solução permanece estável com points = 50; difforder = 2
.
difforder = 4
também pode ser usado se desejar.
Adendo: Reimplementação do Método OP
Depois de olhar mais de perto o código do OP junto com o artigo vinculado no início desta resposta, acho que entendi o que o OP implementou. A seguir está minha implementação para o mesmo método, que pode ser um pouco mais legível:
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]
Mais alguma explicação: neste método, o $\frac{\partial^4}{\partial x^4}$ foi discretizado de forma recursiva, ou seja, $\frac{\partial}{\partial x}$é discretizado primeiro ( C1[[All, All, 1]]
no código do OP, d1
no meu código) e o discretizado$\frac{\partial^4}{\partial x^4}$é calculado usando Dot
. Ainda se sente desconfiado? OK, vamos validar:
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'''"}]

Desde a $\frac{\partial}{\partial x}$, $\frac{\partial^2}{\partial x^2}$ e $\frac{\partial^3}{\partial x^3}$ todos apareceram como intermediários no método, o bcs do problema de OP pode ser imposto modificando a matriz diretamente, por exemplo:
ListLinePlot[{grid, d1withbc . f[grid]}\[Transpose]]

Conforme ilustrado acima, $\frac{\partial f}{\partial x}\Big|_{x=0}=0$ foi imposto.
Uma vez que este código é um DQM de implementação para viga cantilever, então precisamos colocar a condição de contorno correta para tornar este código estável com o número de pontos de grade Np
mudando. Esta é apenas uma pequena modificação, mas funciona para qualquer um Np
, por exemplo
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}]]

Com esta abordagem, temos que considerar Xb[t]
como força externa aplicada ao ponto arbitrário next
como
Xb[t] = 2 G1 Cos[
2 Ω t] (Subscript[x, Np][
t] + α Subscript[x, Np][t]^3) Table[
KroneckerDelta[next, i], {i, 2, Np}];
No caso de next=Np
temos o código acima. A principal razão pela qual o código original produz uma solução instável é a K1M
definição da matriz tirada do artigo Aplicação do método da quadratura diferencial generalizada ao estudo dos fenômenos pull-in de interruptores MEMS, de Hamed Sadeghian, Ghader Rezazadeh e Peter M. Osterberg. Podemos redefini-lo da seguinte maneira
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}];
Podemos comparar esta solução em Xb=0
(pontos vermelhos) com a solução gerada pelo código xzczd com points=10
(linha sólida)

Agora, se colocarmos Np=30
e aplicarmos Xb
ao primeiro ponto como no código acima, obteremos a imagem para cada ponto da grade como segue
Table[Plot[Evaluate[Subscript[x, i][t] /. First@s], {t, 0, tmax},
PlotRange -> All, PlotLabel -> Row[{"i = ", i}]], {i, 1, Np}]

Isso é muito comum responder à força externa. Usando esta matriz K1M = C0 W4
, podemos perceber a ideia principal de Xb
implementação como$x_1(t)$ do seguinte modo
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
e$x_1(t)$diferem por uma constante de cerca de 0,3. Podemos incluir essa constante na condição inicial para$x_1(t)$ mas poderia ser melhor ficar com $x_1(0)=0$como no código acima. Também devemos notar que o algoritmo proposto não é estável para arbitrário Np
, mas podemos torná-lo estável aumentando$\mu$ para o ponto de fronteira $x_1$ como normalmente fazíamos no método das linhas.
{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]

O primeiro canal oscila mais rápido e a amplitude está aumentando exponencialmente. Cada opção de método para metas ou precisão é capaz de computar as oscilações com enorme poder de forma que todos os outros canais cresçam apenas exponencialmente. Em uma forma de intervalo que calcula essa constante, há oscilações.
A otimização é feita com a perspectiva do maior domínio para a solução. Uma vez que todos os canais de solução são dominados por$x_{1}$ isso é o mais importante.
Cortar o domínio permite uma visão informativa:
ReImPlot[%183[[1, 2]], {t, 0, 4.3}, PlotRange -> Full,
ColorFunction -> "DarkRainbow",
PlotLabels -> {"Expressions", "ReIm"}]

A solução de $x_{1}$consistem em uma oscilação mais lenta com uma frequência dependente do tempo que a frequência fica mais rápida com o tempo. Ele oscila abaixo deste envelope mais lento com uma desaceleração lenta com o tempo, mas uma frequência muito mais rápida.
O enredo é impreciso por causa de inferências até ruído no enredo. O ColorFunction
mostra que a oscilação vai até zero. O envelope é assimétrico nas amplitudes em relação ao eixo x.
É uma chance de que a singularidade em 7.12986 e um pouco mais tarde possa ser calculada de forma estável com uma metodologia aprimorada.
As melhores abordagens são
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]]

Entre os dois métodos, há apenas uma pequena diferença, ambos são de alta precisão. Mas as soluções são diferentes. Ambos calculam ruído e erro, no máximo, das oscilações rápidas. Porém, quanto menores as soluções avançam no tempo, mais erros e ruídos se acumulam.
A extrapolação diverge muito mais rápido no momento crítico, 7.12986
mas em subintervalos, as soluções nos outros canais oscilam menos. Uma subdivisão do domínio pode levar a menos oscilação devido à curvatura menos acumulada, supondo que os limites sejam tomados com cuidado. Há uma chance de integrar menos ruído e erro suavizando a oscilação adotando a extrapolação.
Meu problema é que o Método de "Extrapolação" para NDSolve está incompleto nos exemplos. O Mathematica, por outro lado, faz muito mais internamente. Isso também pode ser devido ao alto grau de similaridade entre os dois grupos de métodos apresentados. O cálculo é muito rápido. Existe um ótimo WorkingPrecision
. Isso não pode ser melhorado ainda mais. O comprimento do domínio tem um valor ótimo. Isso me deixa cético.
Eu tenho o conceito de que é apenas um pulso de altura finita e que a curva se acalma em um processo de aniquilação de energia para baixo. Não há nenhum evento catastrófico pela frente. Mas a divergência é muito rápida em muitos pedidos em passos muito pequenos. A computação sempre termina é a mensagem semelhante à mensagem de rigidez de que o tamanho do passo fica muito pequeno. Isso não pode ser superado evitando-se a troca inadequada de rigidez.
A integração adequada de todas as pequenas oscilações pode exigir muito mais tempo e capacidade de computação do que posso oferecer para esta resposta.
A vantagem na primeira parte do domínio computado é bem visualizada por:

As soluções extrapoladas são muito menos oscilantes no subintervalo mais linear. Ele ganha as mesmas oscilações no início e no subintervalo maior que$⁄pi$. O momento de oscilação é muito maior no limite superior do domínio do que com o Stiffnessswitching. Esta é uma comparação da solução selecionada na pergunta.
A avaliação do StepDataPlot
mostra que nestes subintervalos ocorre a mudança de rigidez. No meio, nenhuma rigidez é executada. Isso torna esses cálculos muito mais rápidos do que os da pergunta.
A questão estratégica é se as oscilações em amplitude $-30$são considerados erro / ruído ou fazem parte da solução. Uma vez que o método Runge-Kutta é projetado para zerar o erro e o ruído, a questão é bastante importante. O resultado pode ser transformado na ideia de que a computação em subintervalos é uma otimização da computação no intervalo completo.
O NDSolve já faz esse cancelamento internamente, até certo ponto. Na literatura, esse fenômeno pode ser denominado arco-íris ou caminho para o caos com divergência. Como pode ser verificado no controle de evento programado do tamanho da etapa, essa abordagem não funciona. Ele foi adaptado de uma pergunta em que o NDSolve está operando em uma solução com muitas ramificações. Ele não detectou nenhum ramo.
A subdivisão é provavelmente a melhor, especialmente se a amplitude for absolutamente maior do que $15$. Os valores para os experimentos numéricos são retirados da pergunta. O mais provável é que haja mais interesse.
Para realizar algumas pesquisas sobre o que isso realmente significa, observe a compreensão do método para NDSolve :
Selecione [Flatten [Trace [NDSolve [sistema], TraceInternal -> True]],! FreeQ [#, Método | NDSolve`MethodData] &]
Pergunte a si mesmo: Quais são os métodos da função NDSolve do Wolfram Mathematica?
"Adams" - preditor - corretor Método de Adams com ordens de 1 a 12 "BDF" - Engrenagem fórmulas de diferenciação reversa implícitas com ordens de 1
a 5 "ExplicitRungeKutta" - pares incorporados adaptativos de 2 (1) a 9 (8) métodos de Runge - Kutta " ImplicitRungeKutta "- famílias de Runge implícito de ordem arbitrária - métodos Kutta" SymplecticPartitionedRungeKutta "- Runge intercalado - Métodos Kutta para sistemas hamiltonianos separáveis" MethodOfLines "- método das linhas para solução de PDEs" Extrapolação "- (Gragg -) Método de extrapolação Bulirsch , com possíveis submethods [Bullet] "ExplicitEuler" - método de Euler para a frente [Bullet] "ExplicitMidpoint" - método de regra de ponto médio [Bullet] "ExplicitModifiedMidpoint" - método de regra de ponto médio com suavização de Gragg [Bullet] "LinearlyImplicitEuler" - método de Euler linearmente implícito [Bullet ] "LinearlyImplicitMidpoint" - método de regra de ponto médio linearmente implícito [Bullet] "LinearlyImplicitModifiedMidpoint" - Bader linearmente implícito - método de regra de ponto médio suavizado "DoubleStep" - versão "baby" de "Extrapolation" "LocallyExact" - aproximação numérica para solução simbólica localmente exata "StiffnessSwitching" - permite alternar entre métodos não rígidos e rígidos no meio da
integração "Projection" - invariante - método de preservação "OrthogonalProjection "- método que preserva a ortonormalidade das soluções" IDA "- solucionador de propósito geral para o problema do valor inicial para
sistemas de equações diferenciais - equações algébricas (DAEs)" Tiro "- método de tiro para BVPs" Chasing "- Gelfand - Método de caça Lokutsiyevskii para BVPs" EventLocator "- local do evento para detectar descontinuidades, períodos, etc." FiniteElements "- problemas de elementos finitos
Use Monitoramento e Seleção de Algoritmos :
try[m_] :=
Block[{s = e = 0},
NDSolve[system, StepMonitor :> s++, EvaluationMonitor :> e++,
Method -> m]; {s, e}]
com o Método e opções de real interesse e boas soluções. É triste que este tutorial vá realmente a fundo. O processo de seleção consome muito tempo.
Esta demonstração mostra a metodologia de favor: Plotagem 3D Adaptável na tarefa de plotar uma superfície 3D. Esta é uma demonstração do próprio Stephen Wolfram. E há mais disso. Existe um para plotagem xy: plotagem adaptativa . Este tutorial mostra o método "DoubleStep" para NDSolve . Ele mostra por que o método Runge-Kutta é bem-sucedido para esse problema. Este tutorial leva o leitor de alguma forma ao complexo oculto por trás da Method
opção "Automatic"
que é tão onipresente na solução NDSolve na literatura, a documentação do Mathematica. A boa prática é como obter amostragem adaptativa como na função de gráfico .
O problema é semelhante ao denotado por "para NIntegrate, você deve forçar a avaliação numérica, caso contrário, ele pode empregar algum esquema de quadratura que minimiza o número de pontos de avaliação."
e
"A partir da forma simbólica do integrando NIntegrate pode detectar algumas de suas características, como periodicidade, a fim de minimizar o número de pontos de avaliação. AFAIK aplicará simbólico até que você desligue com Método -> {Automático," Processamento Simbólico "-> Nenhum} (em vez de Automático pode ser a especificação de método explícita) ou usando o método "caixa preta" (_? NumericQ). Ambas as formas não desabilitam o esquema de quadratura. "
Um bom conceito para uma subdivisão é dado em aceleração de gráficos de contorno, amostragem adaptativa para funções de computação lenta em 2d desta comunidade.
O problema fornecido com os dados fornecidos não é rígido, mas fica rígido se a opção de precisão da questão for considerada tão rígida. Como pode ser confirmado estudando a documentação do Mathematica a escolha da recomendação é única WorkingPrecision
.
Descubra como unir várias instâncias da função de interpolação ! O passo importante à frente é que cada período completo deve ser devidamente levado em consideração. Um bom método para subdivisão está documentado no NDSolve EventLocator