# Hyperbolic PDE with spiky behaviour

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).

• Look for PeriodicBoundaryCondition Jul 3, 2023 at 12:23
• I am not quite sure what the suggestion is here. I tried using PeriodicBoundaryCondition and I got worse results. Jul 3, 2023 at 13:00
• 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}]]  Jul 3, 2023 at 13:57
• 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. Jul 3, 2023 at 14:30
• What a pity! Were do the eqs come from? Especially part Exp[2 P[\[Tau], x]]? Jul 3, 2023 at 15:41

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


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

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


• 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. Jul 11, 2023 at 19:21
• 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'' Jul 11, 2023 at 19:23
• @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? Jul 12, 2023 at 3:47
• 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. Jul 12, 2023 at 17:58
• @user12588 Do you mean that picture for A=10 should be like picture for A=8/Pi` shown above? Jul 13, 2023 at 3:49