# Solving wave equation with singular initial conditions

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.

• – anderstood Apr 21 '17 at 18:11

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

• By discretizing with only 76 points, you smoothen the singularity, in a way. For sure it's better, but far from perfect. Do you think it is simply not possible to solve this problem significantly better? – anderstood Apr 21 '17 at 17:12
• @anderstood How do you define "perfect"? If you mean the resulting InterpolatingFunction should involve a non-differentiable point at the cusp, then as far as I can tell it's very hard if not impossible, because 1. currently PDE solver of NDSolve doesn't have a strong enough support for weak solution problem; 2. There's no documented method for the interpolation of weak solution data. (Check this post for more information. ) – xzczd Apr 22 '17 at 8:00
• @xzczd By "far from perfect" I mean "far from exact". The question is about reducing the error. The first version of this answer (when I commented) reduces error, but it is still large: it's easy to see the difference with the exact solution. The edit gives a nice workaround. Your comment also provides relevant information which IMHO would deserve an answer! – anderstood Apr 22 '17 at 23:16

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 : • What is the definition of int? – anderstood Apr 22 '17 at 23:17
• "int" is here for the Integration of v0 from x-ct to x+ct. Here it is just the d'Alembert formula according to initial conditions y0 and v0. It is possible to implement it by several means in mathematica. – CharlelieB Apr 24 '17 at 12:58