6
$\begingroup$

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.

enter image description here

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.

$\endgroup$
1

2 Answers 2

2
$\begingroup$

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

enter image description here

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

enter image description here

Table[Plot[sol[x, t], {x, 0, L}, PlotRange -> All], {t, 0, 1}]
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – anderstood
    Apr 21, 2017 at 17:12
  • 2
    $\begingroup$ @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. ) $\endgroup$
    – xzczd
    Apr 22, 2017 at 8:00
  • $\begingroup$ @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! $\endgroup$
    – anderstood
    Apr 22, 2017 at 23:16
0
$\begingroup$

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 :

Plot of the solution at two different times, showing the propagation of the angular singularity

$\endgroup$
2
  • $\begingroup$ What is the definition of int? $\endgroup$
    – anderstood
    Apr 22, 2017 at 23:17
  • $\begingroup$ "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. $\endgroup$
    – CharlelieB
    Apr 24, 2017 at 12:58

Your Answer

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

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