2
$\begingroup$

I am trying to solve the following systems of hyperbolic PDEs

eqs = {D[
      P[τ, x], {τ, 
       2}] - (Exp[-2 τ] D[P[τ, x], {x, 2}] + 
       Exp[2 P[τ, x]] (D[Q[τ, x], τ]^2 - 
          Exp[-2 τ] D[Q[τ, x], x]^2)) == 0,
   D[Q[τ, x], {τ, 
       2}] - (Exp[-2 τ] D[Q[τ, x], {x, 2}] - 
       2 (D[P[τ, x], τ] D[Q[τ, x], τ] - 
          Exp[-2 τ] D[Q[τ, x], x] D[P[τ, x], x])) == 0};

Subject to the following boundary conditions

bcs = {
\!\(\*SuperscriptBox[\(P\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[0, x] == 10 Cos[x], 
\!\(\*SuperscriptBox[\(Q\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[0, x] == 0, P[0, x] == 0, 
   Q[0, x] == Cos[x], P[τ, 0] == P[τ, 2 π], 
   Q[τ, 0] == Q[τ, 2 π]};

So far, I have tried something like

eqtot = Join[eqs, bcs];
sol = NDSolveValue[
    eqtot, {P[τ, x], Q[τ, x]}, {τ, 0, 6 π}, {x, 0, 
     2 π}, MaxStepSize -> 1/200]; // Timing

However, I am finding poor convergence as a change the MaxStepSize. So, I am searching for ways to improve the code above, and if possible, to find ways to evolve the system to larger times (say $48\pi$ or longer).

$\endgroup$
6
  • $\begingroup$ Look for PeriodicBoundaryCondition $\endgroup$ Jul 3, 2023 at 12:23
  • $\begingroup$ I am not quite sure what the suggestion is here. I tried using PeriodicBoundaryCondition and I got worse results. $\endgroup$
    – user12588
    Jul 3, 2023 at 13:00
  • $\begingroup$ As far as I expirienced the constraint P[\[Tau], 0] == P[\[Tau], 2 \[Pi]] isn't sufficient, one has to use two periodic boundary conditions instead: PeriodicBoundaryCondition[P[\[Tau], x], x == 0, TranslationTransform[{ 2 Pi}]], PeriodicBoundaryCondition[P[\[Tau], x], x == 2 Pi , TranslationTransform[{ -2 Pi}]] $\endgroup$ Jul 3, 2023 at 13:57
  • $\begingroup$ Thank you for your reply. When I implement what you suggest, the code develops a singularity in finite time $\tau$. I cannot see any improvement. $\endgroup$
    – user12588
    Jul 3, 2023 at 14:30
  • $\begingroup$ What a pity! Were do the eqs come from? Especially part Exp[2 P[\[Tau], x]]? $\endgroup$ Jul 3, 2023 at 15:41

1 Answer 1

3
+25
$\begingroup$

This problem can be solved using "MethodOfLines" with options "DifferenceOrder" -> "Pseudospectral". First we rescale equations using substitutions $\tau \rightarrow 2 \pi L t$, and $x\rightarrow 2\pi L x$, then system end boundary conditions take a form

L = 24;
eqs = {D[
      P[t, x], {t, 2}] - (Exp[-4 Pi L t] D[P[t, x], {x, 2}] + 
       Exp[2 P[t, x]] (D[Q[t, x], t]^2 - 
          Exp[-4 Pi L t] D[Q[t, x], x]^2)) == 0, 
   D[Q[t, x], {t, 2}] - (Exp[-4 Pi L t] D[Q[t, x], {x, 2}] - 
       2 (D[P[t, x], t] D[Q[t, x], t] - 
          Exp[-4 Pi L t] D[Q[t, x], x] D[P[t, x], x])) == 0};
 bcs = {
  Derivative[1, 0][P][0, x]/(2 Pi L) == A Cos[2 Pi L x], 
  Derivative[1, 0][Q][0, x] == 0, P[0, x] == 0, 
     Q[0, x] == Cos[2 Pi L x], P[t, 0] == P[t, 1/L], 
     Q[t, 0] == Q[t, 1/L]};
eqtot = Join[eqs, bcs];

Numerical solution

sol = NDSolveValue[eqtot/.A->5/Pi, {P, Q}, {t, 0, 1}, {x, 0, 1/L}, 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", 
      "MaxPoints" -> Round[340 Pi], "MinPoints" -> Round[40 Pi], 
      "DifferenceOrder" -> "Pseudospectral"}}];

Visualization

Table[Plot3D[
  sol[[i]][t/(2 Pi L), x/(2 Pi L)], {t, 0, 48 Pi}, {x, 0, 2 Pi}, 
  PlotPoints -> 100, ColorFunction -> Hue, Mesh -> None, 
  MaxRecursion -> 2, AxesLabel -> {"\[Tau]", "x"}, 
  PlotRange -> All], {i, 2}]

Figure 1

Then we can play with parameter A to get stably result. For instance at A=8/Pi we have in 0.15 s this nice picture Figure 2

Update 1. This is relatively stably result for A=10

eqs = {D[
      P[t, x], {t, 2}] - (Exp[-4 Pi L t] D[P[t, x], {x, 2}] + 
       Exp[2 P[t, x]] (D[Q[t, x], t]^2 - 
          Exp[-4 Pi L t] D[Q[t, x], x]^2)) == 0, 
   D[Q[t, x], {t, 2}] - (Exp[-4 Pi L t] D[Q[t, x], {x, 2}] - 
       2 (D[P[t, x], t] D[Q[t, x], t] - 
          Exp[-4 Pi L t] D[Q[t, x], x] D[P[t, x], x])) == 0};
L = 24; bcs = {
  Derivative[1, 0][P][0, x]/(2 Pi L) == A Cos[2 Pi L x], 
  Derivative[1, 0][Q][0, x] == 0, P[0, x] == 0, 
     Q[0, x] == Cos[2 Pi L x], P[t, 0] == P[t, 1/L], 
     Q[t, 0] == Q[t, 1/L]};
eqtot = Join[eqs, bcs];
sol = NDSolveValue[eqtot /. A -> 10, {P, Q}, {t, 0, 1}, {x, 0, 1/L}, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MaxPoints" -> Round[120 Pi], "MinPoints" -> 81, 
        "DifferenceOrder" -> 4}}]; // AbsoluteTiming

Figure 3

$\endgroup$
6
  • $\begingroup$ I don't think this is quite right. I think you forgot to change the BCs when you did the replacement $\tau\to2\pi L t$ and $x\to2\pi L x$. In particular, the initial condition for the derivative of $P$ gets a factor of $L$, which changes the behaviour dramatically. In particular, I again see no convergence. $\endgroup$
    – user12588
    Jul 11, 2023 at 19:21
  • $\begingroup$ The best I was able to do was to do the following ``n = 7000; eqtot = Join[eqs, bcs]; sol1 = NDSolveValue[ eqtot, {P[[Tau], x], Q[[Tau], x]}, {[Tau], 0, 12 [Pi]}, {x, 0, 2 [Pi]}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "TemporalVariable" -> [Tau], "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> n, "MaxPoints" -> n}}}]; // Timing'' $\endgroup$
    – user12588
    Jul 11, 2023 at 19:23
  • $\begingroup$ @user12588 You are right. Code been updated. Maybe we can play with parameter A in Derivative[1, 0][P][0, x]/(2 Pi L) == A Cos[2 Pi L x] to get stable result. Is your current parameter A=10 related to some physical reason or just to your arbitrary choice? $\endgroup$ Jul 12, 2023 at 3:47
  • $\begingroup$ Alas, I don't think the plot generated by your latest attempt is very stable. The one I posted previously does a better job, I think. My feeling is that something like Discontinuous Galerkin Method (DGM) would be useful, possibly with some sort of dynamical mesh refinement as the spikes develop. There are a number of packages (for Python or Fortran users) that do this, such as `Clawpack' (clawpack.org). These are old enough that I was hoping that MMA would have these readily available. $\endgroup$
    – user12588
    Jul 12, 2023 at 17:58
  • $\begingroup$ @user12588 Do you mean that picture for A=10 should be like picture for A=8/Pi shown above? $\endgroup$ Jul 13, 2023 at 3:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.