Solving wave equation with singular initial conditions

Question

I am trying to compute the solution to the wave equation using NDSolveValue and singular initial displacement (continuous, but not differentiable in a point). I need to propagate as accurately as possible the singularity of the initial condition through time. However the singularity is "smoothened" during the integration, see below.

The red curve is my initial condition in displacement. The purple one is the output at a different time of the computation. As you can see, the purple curve is perfectly smooth, while it shouldn't.

Here you can see a piece of the script I use to get these solutions :

(*initialisation of parameters and initial conditions*)
Sigma0 = 0.002; g0 = 0.001;  L = 1; ρ = 1; Ei = 1; c =Sqrt[Ei/ρ]; epsi = 10^-7.;
qu[x_] := 
Piecewise[{{(-1 Sigma0/Ei) x, 0 <= x <= L/2}, {(-1 Sigma0/Ei) L/2, 
 L/2 <= x <= L}}];
qv[x_] := 
  Piecewise[{{0, 0 <= x <= L/2}, {(c Sigma0/Ei), L/2 <= x <= L}}];

(*Building of a periodic Solution according to the wave equation 
  in Dirichlet-Neumann Boundary condition*)

u0[x_] :=  qu[x];
v0[x_] :=  qv[x];
uif = NDSolveValue[{D[u[x, t], {t, 2}] - D[u[x, t], {x, 2}] == 0, 
    u[x, 0] == u0[x], (D[u[x, t], t] /. t -> 0) == 
     v0[x], (D[u[x, t], x] /. x -> L) == 0, u[0, t] == 0}, 
   u, {x, 0, L}, {t, -10, 10}];
tc = t /. FindRoot[uif[L, t] - g0, {t, 0.1}];

(*Plot of the solution at a time tc in purple and the initial condition*)
Show[Plot[uif[x, 0], {x, 0, L}, PlotPoints -> 200, PlotStyle -> Red], 
 Plot[uif[x, tc], {x, 0, L}, PlotPoints -> 200, PlotStyle -> Purple], 
 PlotRange -> {{0, L}, All}, AxesLabel -> {"Time", "Displacement"}]

I get those warnings during all the computation :

NDSolveValue::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x.

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 185.53982144076966 at t = 20. in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

NDSolveValue::ndstf: At t == -2.61883, system appears to be stiff. Methods Automatic, BDF, or StiffnessSwitching may be more appropriate.

Answer 1

NDSolve solution

The problem is that NDSolve can not handle non-smooth initial condition very well.

An alternative, is to change the method of solution from default to MethodOfLines

uif = NDSolveValue[{D[u[x, t], {t, 2}] - D[u[x, t], {x, 2}] == 0, 
     u[x, 0] == u0[x], (D[u[x, t], t] /. t -> 0) == v0[x], (D[u[x, t], x] /. x -> L) == 0, 
     u[0, t] == 0}, u, {x, 0, L}, {t, -10, 10}, 
     Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", 
     "MinPoints" -> 5*15 + 1, "MaxPoints" -> 5*15 + 1,"DifferenceOrder" -> Automatic}}];

tc = t /. FindRoot[uif[L, t] - g0, {t, 0.1}];

Show[Plot[uif[x, 0], {x, 0, L}, PlotPoints -> 200, PlotStyle -> Red], 
 Plot[uif[x, tc], {x, 0, L}, PlotPoints -> 200, PlotStyle -> Purple], 
 PlotRange -> {{0, L}, All}, AxesLabel -> {"Time", "Displacement"}]

DSolve solution

An analytical solution can be obtained using DSolveValue

sol = DSolveValue[{D[u[x, t], {t, 2}] - D[u[x, t], {x, 2}] == 0, 
   u[x, 0] == u0[x], (D[u[x, t], t] /. t -> 0) == v0[x]}, u, {x, 0, L}, {t, -10, 10}]

Table[Plot[sol[x, t], {x, 0, L}, PlotRange -> All], {t, 0, 1}]

Answer 2

I have found a way to avoid the solver. For the wave equation, it is possible to propagate effectively the singularities with the d'Alembert's solution. It works with this case but it cannot be implemented with other PDE. I used my initial conditions to create a solution according to this formula :

y[x,t]=1/2y0[x-ct]+1/2y0[x+ct]+1/(2c)int[v0[s],{(x-ct),(x+ct)}]

Here is my result :