I have two extra points I would like to add to Wolfgang's answer. A formulation
of the CFL condition that I find more useful than the classic formula is this:
A necessary condition for the stability of a numerical scheme is that the
numerical domain of dependence bounds the physical domain of dependence.
This is exactly what good old
$$
\dfrac{\Delta t}{c \Delta x} \leq 1
$$
says for simple schemes. The quote I wrote above has a bit more nuance, though,
that can be used to answer the questions about implicit methods and higher order
algorithms (my two points).
My points below are for the first order wave equation just to make things easy:
I hope that everyone can see how the same ideas apply to the second order
equation.
Higher order algorithms
It is true that for most time stepping methods increasing the order imposes
stricter stability constraints but a notable exception to this is explicit Runge
Kutta methods (which include some part of the imaginary axis in their stability
regions; that is a separate stability concern). Consider the first order wave
equation $$ u_t + c u_x = 0 $$ with periodic boundary conditions and second
order finite differences $$ u_x(x_j) = \dfrac{u(x_j + \Delta x) - u(x_j - \Delta
x)}{2 \Delta x} + O(\Delta x^2) $$
(to make things as easy as possible) in space. One can write the semidiscrete
system as
$$
\dfrac{d\vec{U}}{dt} = A \vec{U}.
$$
If we use RK3 then we have
$$
\begin{align}
k_1 &= A \vec{U}\,^n \\
k_2 &= A \left(\vec{U}\,^n + \dfrac{\Delta t}{2} k_1\right) \\
k_3 &= A \left(\vec{U}\,^n + \dfrac{3 \Delta t}{4} k_2\right) \\
\vec{U}\,^{n + 1} &= \vec{U}\,^n + \dfrac{\Delta t}{9} (2 k_1 + 3 k_2 + 4
k_3) \\
\end{align}
$$
so, by definition of (A):
- the $j$th entry of $k_1$ uses $U_{j - 1}^n$ and $U_{j + 1}^n$
- the $j$th entry of $k_2$ uses $U_{j - 2}^n$, $U_{j - 1}^n$, $U_j^n$, $U_{j + 1}^n$, and $U_{j + 2}^n$
- the $j$th entry of $k_3$ uses $U_{j - 3}^n$ through $U_{j + 3}^n$
since for each stage we use the same space operator on the current solution
which effectively widens the stencil. Hence doing one step with this method has
the CFL constraint
$$
\dfrac{c \Delta t}{\Delta x} \leq 3
$$
since calculating $U_j^{n + 1}$ requires data from three grid points to the left
and three grid points to the right. The same argument holds for higher order
methods: a method with more stages tends to have a wider stencil and therefore a
larger Courant number.
Warning: due to the spatial operator we used this scheme is not stable, even
if it satisfies the CFL constraint. The correct Courant number to use (which
comes from the stability region of RK3) is (roughly)
$$
\dfrac{c \Delta t}{\Delta x} \leq 1.75.
$$
Implicit methods
The answer to your question about ''well-defined values for $C_{max}$ when using
implicit methods'' is that $C_{max} = \infty$ by the same plain-text version of
the CFL condition. To see why this is true consider the same ODE system with
backward Euler:
$$
\vec{U}\,^{n + 1} = \vec{U}\,^{n} + \Delta t A \vec{U}\,^{n + 1}
\Rightarrow
(I - \Delta t A) \vec{U}\,^{n + 1} = \vec{U}\,^{n}
$$
so solving for the solution at $t^{n + 1}$ requires knowledge of all values at
$t^n$, meaning that no matter what time step we use we are pulling in the
entire physical domain of dependence. Implicit methods for wave equations are an
unusual choice for a variety of (unrelated to the original question) reasons but
they do satisfy the CFL condition for any choice of $\Delta t > 0$.
Edit:
@Mathews24 asked some more closely related questions whose answers won't
really fit in a comment so I am appending them to my answer here. I hope that
that is an acceptable way to answer questions (please correct me if there is a
better place for this).
Why is the centered difference scheme unstable past a Courant number of $1.75$?
As mentioned previously the CFL condition is a necessary, but not sufficient,
condition for the stability of a numerical scheme. The other usual stability
condition comes from the ODE solution method (which depends on the eigenvalues
of $A$, which in turn depends on the spatial discretization: this is why I say
that the choice of spatial discretization impacts the stability of the time
discretization). Ignoring the physical meaning of things for a moment, we say
that an ODE solver for a scalar problem
$$
y' = \lambda y + f
$$
is stable if the numerical approximation of the unforced problem, $y' =
\lambda y$, remains bounded for all time. I like to think of this
thermodynamically: if one considers an unforced wave in a closed system then
that wave will not grow in energy over time (but we allow it to shrink for the
purposes of our numerical scheme, since that is better than exponential
growth). An easy example of this is forward Euler: one can check that if we do
$$
w^n \approx y(t^n), \,w^{n + 1} = w^n + \lambda \Delta t w^n
\Rightarrow w^{n + 1} = (1 + \lambda \Delta t) w^n
$$
that the scheme is stable when (for complex $\lambda$) $|\lambda \Delta t + 1|
\leq 1$. The value on the left is called the amplification factor: the
solution will grow exponentially if it is higher than one. The same idea holds
for all other ODE schemes (in particular Runge-Kutta methods): one can calculate
a stability region for the three stage scheme I presented above with a
restriction on $\lambda \Delta t$ that guarantees that the solution does not
grow. The stability plot is the blue curve under ''Bogacki-Shampine-4-2-3'' at
http://runge.math.smu.edu/arkode_dev/doc/guide/build/html/Butcher.html
if anyone wants to see it.
The rule of thumb that people use (that works for all non-pathological cases;
there is a nice book by Trefethen and Reddy about how to do this more generally)
to examine stability of a system of ODEs $y' = A y$ is to do 1D stability
analysis with all of the eigenvalues of $A$. This is rigorous if the matrix is
diagonalizable and usually works even if it is not. If one looks at the
stability region of the RK3 scheme then one will see that it stays (barely!) to
the right of the imaginary axis and intersects at about $1.75$. For the matrix I
gave above the eigenvalues are purely imaginary and are bounded (in modulus) by
$$
|\text{eig}(A)| \leq \dfrac{c}{\Delta x}
$$
which is exactly why things blow up for a Courant number greater than $1.75$.
Why are implicit methods unusual for wave equations?
Running this problem with RK3 near the stability limit is nice because the
energy in the system is very nearly conserved: things decay very slightly over
time since all of the amplification factors are very close to $1$ (the product
$\lambda_i \Delta t$ is always near the boundary of the stable region). If we
used an implicit method then we won't have to worry about the CFL constraint
(or, for the right choice of implicit method, the amplification factors) but one
can check that as a result of this the amplification factors will be
significantly less than one, leading to a lot of extra dissipation in the
scheme: this is usually not what you want with waves (usually some sort of
conservation is desired, even if it is only approximate). I recommend
calculating the amplification factors of backward Euler along the imaginary axis
to convince oneself of why this is a concern.
