Consider the damped wave equation with Dirichlet boundary conditions
$\frac{\partial^2}{\partial t^2}+2b(x)\frac{\partial}{\partial t}=\frac{\partial^2 u}{\partial x^2}$, $0<x<\pi $, $0\leq t<\infty$, $b(x)>0 \forall x$, $ u(x,0)=u_{0}(x)$, $\frac{u(x,0)}{\partial t}=u_{1}(x)$, $u(0,t)=u(\pi,t)$.
Smoothness is assumed for $b(x)$, $u_{0}$, $u_{1}$ and $u$.
One of the questions that it is worth asking here is the decay rate of the total energy of the system. For constant b, it is trivial to calculate it via separation of variables, however, for varying $b(x)$ it is not straightforward.
The simplest method to try to do that is via a semi-discretization, where we divide the space in N nodes, each node satisfying the equation
$$\ddot{u}_{j}+2b_{j}\dot{u}_{j}=\frac{-u_{j+1}-u_{j-1}+2u_{j}}{(\Delta x)^2} $$
where we used the second-order approximation for the second derivative in order to get it, and divided the damping in $N$ parts.
These equations form a system of PDEs whose Jacobian is
$ \begin{matrix} 0 & 1 \\ -L & -B \end{matrix}$, $B_{ij}=\delta_{ij}b_{i}$,
and $L$ is the Laplacian matrix of an array plus 1 in the elements $(i,j)=(1,1)$, $(i,j)=(N,N)$.
The eigenvalue with the highest real part would correspond then to the decay rate.
One would think naively that for high enough $N$ that this approximation leads to the correct decay rate. However, it is known that this is not the case, as one can check in the paper 'Designing for optimal energy absorption, III: numerical minimization of the spectral abscissa', in the section 3.
This is mainly do to the fact that even for high $N$, the highest eigenvalues of $-L$ will NOT converge to those of $\frac{d^2}{dx^2}$, leading to errors in the eigenvalues of the whole system. This results seem to highlight that using such approximation should not lead to a convergence of the solution, i.e., the difference between $u_{approx}$ and $u$ is a function of $\Delta x$ that goes to zero as $\Delta x$ goes to 0.
However, at the same time, we know that a finite element analysis of the system (i.e. discretizing both the time and the position) leads to a numerical solution that converges to the actual solution when the discretization over time and position is small.
Discretizing over time and position is equivalent from discretizating first in relation to the position, and then in relation to the time. Therefore, I believe that if one method converges, the other should also converge.
What am I missing here?