# Why does the FDM give a correct solution to a PDE with a discontinuous initial condition?

I was solving the dimensionless wave equation:

$$u_{xx} = u_{tt} \tag 1$$

with the initial conditions:

$$u(x,0) = 0 \tag 2$$

$$u_t(x=0,0) = v_0 \tag 2$$ $$u_t(x>0,0) = 0 \tag 3$$

and the boundary conditions:

$$u(1,t) = 0 \tag 4$$ $$mu_{tt}(0, t) + u_x(0, t) = 0 \Leftrightarrow mu_{xx}(0, t) + u_x(0, t) = 0 \tag 5$$

This represents the mathematical model of an elastic rod subjected to a rigid body impact. The constant $$m$$ is the mass ratio of the rigid body and elastic rod, and $$v_0$$ is the velocity of the rigid body at the initial time of impact. The function $$u$$ is the longitudinal displacement.

I managed to solve this model up to $$t_{max} = 2$$. The solution for $$t \in [0, 1]$$ is:

$$u(x,t) = mv_0 \Bigg( e^{ {t-x}\over m } - 1 \Bigg) \ \mathrm{if} \ t \geq x \tag 6$$ $$u(x,t) = 0 \ \mathrm{if} \ t \leq x \tag 7$$

For $$t \in [1, 2]$$, the solution is:

$$u(x,t) = mv_0 \Bigg( e^{ {t-x}\over m } - e^{ {t+x-2}\over m } \Bigg) \ \mathrm{if} \ t \geq x - 2 \tag 8$$ $$u(x,t) = mv_0 \Bigg( e^{ {t-x}\over m } - 1 \Bigg) \ \mathrm{if} \ t \leq x - 2 \tag 9$$

The problem is that the solution is piecewise-defined. The expression changes each time the wave hits a boundary. So using an analytical solution to analyze the wave propagation is not very practical if the impact duration is long enough to have a lot of wave reflections. So I decided to simulate this process using the finite difference method.

To my surprise, the numerical solution and the analytical solution match perfectly, as shown in the graphs below:

Note that the ratio of the time step and spatial step is not equal to one (the so-called magic CFL number). It is equal to $$0.67$$.

Here is what I don't understand. I have a discontinuous initial condition $$(2)$$, $$(3)$$. The group velocity and phase velocity are not equal. I was expecting a significant numerical dispersion at the wavefront. However, that did not happen. My question is, why?

In the past, I was solving a similar problem and had an issue with significant numerical dispersion. I even made a post about it here. Why is it that this time the numerical dispersion is not significant enough to be visible when plotting the solution? Could this be only because my solution is continuous even though it is not smooth (so it's not the smoothness of the solution, but the continuity of it that affects numerical dispersion)?

Note: I will be making edits suggested by comments.

EDIT 1:

I used the central difference approximation of second-order accuracy:

$$u_{tt} = {U_{i+1,j} - 2U_{i,j} + U_{i-1,j}\over \Delta t^2} \tag {10}$$ $$u_{xx} = {U_{i,j+1} - 2U_{i,j} + U_{i,j-1} \over \Delta x^2} \tag {11}$$

EDIT 2:

I calculated analytically and numerically the partial time derivative of displacement. The numerical partial derivative was calculated using the central difference approximation. I came across something that might give a clue to finding the answer to my question. Namely, if the CFL number is equal to one, there is a good match between the solutions (as shown in the two pictures below):

However, after changing the CFL number to $$0.67$$, the displacement wave function still did not deviate noticeably, but the time derivative of it did. Especially at the wavefront (as shown in the picture below). The oscillations are taking place:

EDIT 3:

Because the graphs seem more like straight lines than exponential functions (which is visually misleading), I added a code snippet showing how I calculated the analytical solution in Matlab. Indeed, the graphs are exponential functions, not linear ones.

%% Analytical solution

u   = zeros(nx, not);
u_t = zeros(nx, not);

u_t(1, 1) = v_0;

if t_max == 2

for i = 1 : (nt + 1) / 2

for j = 1 : nx

if t(i) > x(j)

u(j, i)   = v_0 * m * (exp( (t(i) - x(j)) / m ) - 1);

u_t(j, i) = v_0 * exp( (t(i) - x(j)) / m );

end
end
end

for i = ( (nt + 1) / 2) + 1 : not

for j = 1 : nx

if t(i) > (2 - x(j))

u(j, i)   = v_0 * m * (                        ...
exp( ( t(i) - (x(j)    ) ) / m ) - ...
exp( ( t(i) + (x(j) - 2) ) / m )   ...
);

u_t(j, i) = v_0 * (                            ...
exp( ( t(i) - (x(j)    ) ) / m ) - ...
exp( ( t(i) + (x(j) - 2) ) / m )   ...
);

else

u(j, i)   = v_0 * m * (exp( (t(i) - x(j)) / m ) - 1);

u_t(j, i) = v_0 * exp( (t(i) - x(j)) / m );

end
end
end
end


Also, I posted a code snippet below showing how I use finite differences to calculate the wave functions.

%% Numerical functions

U   = zeros(nx, nt);
U_t = zeros(nx, nt);

%% Initial condition implementation

U(1, 2)   = dt * v_0;
U_t(1, 1) = v_0;

%% Finite difference loop

j = 2 : nx - 1;

for i = 2 : nt - 1

g = ( 4 * m * U(1, i) - (2 * m + dx) * U(1 + 1, i) ) / (2 * m - dx); % Ghost node

U(1, i + 1) = -U(1, i - 1) + 2 * U(1, i) + (dt^2 / dx^2) * ( U(1 + 1, i) - 2 * U(1, i) + g );

U(j, i + 1) = -U(j, i - 1) + 2 * U(j, i) + (dt^2 / dx^2) * ( U(j + 1, i) - 2 * U(j, i) + U(j - 1, i) );

U_t(:, i) = ( U(:, i + 1) - U(:, i - 1) ) / (2 * dt);
end

%% Calculate time derivative at last time step

UU = zeros(nx, 1);

i = nt;

g = ( 4 * m * U(1, i) - (2 * m + dx) * U(1 + 1, i) ) / (2 * m - dx);

UU(1) = -U(1, i - 1) + 2 * U(1, i) + (dt^2 / dx^2) * ( U(1 + 1, i) - 2 * U(1, i) + g );

UU(j) = -U(j, i - 1) + 2 * U(j, i) + (dt^2 / dx^2) * ( U(j + 1, i) - 2 * U(j, i) + U(j - 1, i) );

U_t(:, i) = ( UU - U(:, i - 1) ) / (2 * dt);


EDIT 4: Below I added the error graphs suggested in the comment section. I calculated the error as the analytical solution minus the numerical solution ($$Err = u - U$$). In the titles are the time-space ratio, as well as the time moments. When the time-space ratio is one, nu dispersion is noticed, while in the other case, the dispersion is noticeable but of small magnitude.

• You haven't said what numerical scheme you used and that could make a huge difference. To address at least one part of your question, there's a big difference between having a 0th-order and 1st-order time derivative being discontinuous. Commented Mar 10 at 21:36
• I edited my answer. It was the 5-point central difference stencil. Also, I agree with you on the derivative part. However, if $u_t$ is discontinuous, then also $u_{tt}$ and $u_{xx}$ should be discontinuous. I know I am integrating twice in time, but I am also differentiating twice in space which I am not sure why gives a good result. I don't know why differentiating twice in space doesn't produce a problem. Commented Mar 10 at 23:41
• Having made more than my fair share of mistakes, could you check if you are indeed plotting the numerical and analytical solution? Or perhaps, are you plotting the analytical solution twice?
– NNN
Commented Mar 12 at 11:20
• Did you check convergence, i.e., the error in infinity norm when increasing spatial resolution / drecrease time step? Does that match your expectations? Commented Mar 12 at 14:57
• @NikolaRistic I think you should visualize the derivative of your solution. In particular, it is consistent with your claim of using discontinuous initial conditions. Otherwise you would compare apples with oranges. Commented Mar 13 at 9:03

Looking at the plots of the errors, it seems that there is not much difference between the oscillations that you have in this case with respect to those that you had in your previous post (Can this finite difference dispersion be eliminated somehow?). In your previous post the peak of the oscillations was approximately as high as $$1/3$$ of the jump that you have in your solution. Here I see approximately the same.
Let's pick the plots at $$t=0.4$$ for the case where $$\frac{\Delta t}{\Delta x } \neq 1$$: the oscillation peaks at $$0.5$$E-3 and the $$\Delta u$$ is approximately $$0.016$$ in the vicinity of the peak (you have 12 points ranging from $$u=0.2$$ to $$u=0$$). The ratio $$0.005 / 0.016$$ is approximately $$1/3$$.