Discrete energy conservation
In the absence of additional factors (dissipation, external forcing), the total energy should be conserved.
The discrete equations discussed in the temporal discretization are designed to satisfy this property, which is confirmed here.
The simplified discrete equations are
\[\begin{split}\begin{pmatrix}
1 & - \tau \\
C_k \tau & 1
\end{pmatrix}
\begin{pmatrix}
\kpos{k}{n + 1} \\
\kvel{k}{n + 1}
\end{pmatrix}
=
\begin{pmatrix}
1 & \tau \\
- C_k \tau & 1
\end{pmatrix}
\begin{pmatrix}
\kpos{k}{n } \\
\kvel{k}{n }
\end{pmatrix}.\end{split}\]
We define the discrete total kinetic and potential energies as
\[ \begin{align}\begin{aligned}K
&
=
\mysum{k}{0}{N - 1}
L
\frac{1}{2}
\rho
\left( \kvel{k}{} \right)^2,\\U
&
=
\mysum{k}{0}{N - 1}
L
\frac{1}{2}
\rho
C_k
\left( \kpos{k}{} \right)^2,\end{aligned}\end{align} \]
which are identical to the theoretical relations but truncated at \(k = N\).
Their differentiations (in time) yield
\[ \begin{align}\begin{aligned}\dif{K}{t}
&
=
\mysum{k}{0}{N - 1}
L
\frac{1}{2}
\rho
\frac{
\left( \kvel{k}{n + 1} \right)^2
-
\left( \kvel{k}{n } \right)^2
}{\Delta t}\\&
=
\mysum{k}{0}{N - 1}
L
\rho
\frac{
\kvel{k}{n + 1}
-
\kvel{k}{n }
}{\Delta t}
\frac{
\kvel{k}{n + 1}
+
\kvel{k}{n }
}{2}\\&
=
\mysum{k}{0}{N - 1}
L
\rho
\dif{\kvel{k}{}}{t}
\ave{\kvel{k}{}},\\\dif{U}{t}
&
=
\mysum{k}{0}{N - 1}
L
\frac{1}{2}
\rho
C_k
\frac{
\left( \kpos{k}{n + 1} \right)^2
-
\left( \kpos{k}{n } \right)^2
}{\Delta t}\\&
=
\mysum{k}{0}{N - 1}
L
\rho
C_k
\dif{\kpos{k}{}}{t}
\ave{\kpos{k}{}}\\&
=
\mysum{k}{0}{N - 1}
L
\rho
C_k
\ave{\kvel{k}{}}
\,
\ave{\kpos{k}{}},\end{aligned}\end{align} \]
and, as a consequence:
\[ \begin{align}\begin{aligned}\dif{E}{t}
&
=
\dif{K}{t}
+
\dif{U}{t}\\&
=
\mysum{k}{0}{N - 1}
L
\rho
\ave{\kvel{k}{}}
\left(
\dif{\kvel{k}{}}{t}
+
C_k
\ave{\kpos{k}{}}
\right).\end{aligned}\end{align} \]
Since we design the discrete equation of motion such that the term inside the parentheses is zero, we see that our scheme satisfies the discrete energy conservation:
\[\dif{E}{t}
=
0.\]
External forcing
The simplest approach to discretize the term is to adopt the Euler-forward scheme, yielding
\[\myint{
t^n
}{
t^{n+1}
}{
\intfactor{t - t^n}
A_k
}{
t
}
\approx
A_k^{n}
\Delta t.\]
Although external forcing can be arbitrarily prescribed, we limit our focus to a sound source, i.e., a specific profile in physical space is given, whose amplitude only changes in time:
\[\myint{
t^n
}{
t^{n+1}
}{
\intfactor{t - t^n}
A_k
}{
t
}
=
A_k
\myint{
t^n
}{
t^{n+1}
}{
\intfactor{t - t^n}
\sin \left( \omega t \right)
}{
t
}.\]
Due to
\[ \begin{align}\begin{aligned}&
\left(
1
+
\frac{N_k^2}{\omega^2}
\right)
\myint{
t^n
}{
t^{n+1}
}{
\intfactor{t - t^n}
\sin \left( \omega t \right)
}{
t
}\\=
&
\left[
-
\frac{1}{\omega}
\intfactor{t - t^n}
\cos \left( \omega t \right)
+
\frac{N_k}{\omega^2}
\intfactor{t - t^n}
\sin \left( \omega t \right)
\right]_{t^n}^{t^{n + 1}}\\=
&
\intfactor{\Delta t}
\left\{
-
\frac{1}{\omega}
\cos \left( \omega t^{n + 1} \right)
+
\frac{N_k}{\omega^2}
\sin \left( \omega t^{n + 1} \right)
\right\}
-
\left\{
-
\frac{1}{\omega}
\cos \left( \omega t^{n } \right)
+
\frac{N_k}{\omega^2}
\sin \left( \omega t^{n } \right)
\right\},\end{aligned}\end{align} \]
we are able to evaluate it analytically.