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.
