O método de quadratura diferencial falha no PDE de 4ª ordem com não linear bc conforme a grade fica mais densa

Dec 26 2020

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, NDSolvefunciona 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 ( Npno código abaixo) aumenta. Se tmaxfor 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 - 1ODEs 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 = 6e 8sã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 NDSolvepara lidar com o sistema rígido, mas ainda não está funcionando.

Pensei que, se usar NDSolveem 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 NDSolvesimulaçã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

9 xzczd Jan 03 2021 at 11:40

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 Pseudospectralpode 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 NDSolvediretamente 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 pointsaumenta 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, d1no 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.

7 AlexTrounev Dec 28 2020 at 06:04

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 Npmudando. 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 nextcomo

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=Nptemos o código acima. A principal razão pela qual o código original produz uma solução instável é a K1Mdefiniçã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=30e aplicarmos Xbao 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 Xbimplementaçã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}]

Finalmente podemos verificar isso Xbe$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]}  

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]

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 ColorFunctionmostra 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.12986mas 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 StepDataPlotmostra 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 Methodopçã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