For the second-order hyperbolic PDE on $\Omega = (0,1)$:
$$ u_{tt}(x,t) = c^2 u_{xx}(x,t)\; , \quad \begin{cases}
u(0,t) = L(t) \\
u(1,t) = R(t) \\
u(x,0) = f(x) \\
u_t(x,0) = g(x)
\end{cases} $$
I would like to convert it into a first-order system and apply Lax-Friedrichs scheme. In particular, let
$$ r(x,t) = c u_x(x,t) \quad \text{and} \quad s(x,t) = u_t(x,t)$$
so that we have
$$ \begin{cases}
r_t(x,t) = cs_x(x,t) \\
s_t(x,t) = cr_x(x,t)
\end{cases} $$
By applying the Lax-Friedrichs scheme with $\{x_0, x_1, \cdots, x_J\}$ and $\{t^0, t^1, \cdots t^N\}$, we have
$$ \begin{cases}
r^{n+1}_j = \frac{1}{2}\left( r^{n}_{j+1} + r^{n}_{j-1} \right) + \frac{c\Delta t}{2\Delta x} \left( s^{n}_{j+1} – s^{n}_{j-1} \right) \\
s^{n+1}_j = \frac{1}{2}\left( s^{n}_{j+1} + s^{n}_{j-1} \right) + \frac{c\Delta t}{2\Delta x} \left( r^{n}_{j+1} – r^{n}_{j-1} \right) \\
u^{n+1}_j = u^n_j + \Delta t \left(s^n_j\right)
\end{cases} $$
However, it looks like we need boundary and initial conditions for both $r$ and $s$ to solve this system. I know we have
$$ s^0_j = g(x_j) $$
How should I construct the remaining conditions for $r^0_j$, $r^n_0$, $r^n_J$, $s^n_0$ and $s^n_J$?
Best Answer
The system of interest is
$$ \boldsymbol w_t + \begin{pmatrix} 0 & -c \\ -c & 0 \end{pmatrix} \boldsymbol w_x = \boldsymbol 0 \tag{1} \label{1}$$
with initial conditions
$$r(x,0) = c u_x(x, 0) = c \partial_x f(x), \quad s(x,0) = u_t(x, 0) = g(x). $$
The prescription of Dirichlet boundary conditions, i.e., $$u(0, t) = L(t), \quad u(1, t) = R(t) $$ is often not well-posed for hyperbolic systems.
The characteristic curves of the linearized system $\eqref{1}$ are $$ x_\pm = x_0 \pm c t.$$
Consequently, the initial data is transported with speed $c$ to the left and right, respectively. If the functions $L(t), R(t)$ are not constructed in a way such that the problem is well-posed, i.e., that $L(t), R(t)$ correspond to the initial data reaching the boundaries, there exists no solution and you simulation will produce meaningless results.
For general initial data (i.e., except for constant functions) I do not see how you would construct such $L(t), R(t)$, as there are multiple characteristic curves originating from different $x_0$ which reach the boundaries $0, 1$ at each time $t$.
If $r(x, 0) = c u_x = c f_x $ and $s(x, 0) = u_t(x, 0) = g(x)$ are both constant in $x$, however, you can use d'Alemberts Formula to compute the solution of the wave equation as
$$u(x, t) = g + ht = L(t) = R(t) $$ which is independent of space (as expected for initial data that does not depend on space) but still depends on time.
Something that should work, in contrast, are outflow boundary conditions $$r_0 = r_1, \: r_N = r_{N-1}, \quad s_0 = s_1, \: s_N = s_{N-1}$$ which are essentially Neumann boundary conditions that do not clash with the characteristics.