The Wave Equation as a Linear System
In our previous class, we completed our exploration of higher-order explicit time-integration techniques, and moved from a simple spring system to a more complex 1-D Wave Equation. Using the fourth-order Runge-Kutta (RK4) integration scheme, we adapted our spring solver to the wave equation, which seemed stable at first, but as soon as we introduced external forces, the system showed immediate instabilities.What follows in this post is a step-by-step progression, starting with the original PDE of the 1D Wave Equation, and eventually producing a linear system, represented by the equation \(A x = b\), where \(A\) is a square, mostly sparse matrix of size \(N\times N\), where \(N\) is the number of unknown values in the discretized simulation grid (either heights, or accelerations). This post is a lot of math, but all of it is manipulation of simple terms and equations in a linear algebra setting. I've broken things down to the smallest steps because this is such a core concept to how fluid simulations and grid simulations are performed.
The problem in our previous solution is related to how we calculated our fluid vertical acceleration from the fluid's height. To start with, the 1D Wave Equation is written as:
\[
\begin{equation}
\frac{\partial^2 h}{\partial t^2} = c^2 \frac{\partial^2 h}{\partial x^2}
\label{eq:Wave1D}
\end{equation}
\]
Where \(h\) represents height at a given location along an axis, and \(c^2\) represents the square of the wave propagation speed. I think this is my favorite equation, because of the implied equivalency of spatial and temporal dimensions. It says that, as far as waves are concerned, space and time are interchangeable (assuming you choose your units correctly). Anyway, let's examine how we discretized this solution for a given point in our height array, \(h_i\) and its acceleration \(a_i\).
\[
\begin{equation}
a_i = c^2 \frac{h_{i+1} -\ 2 h_i + \ h_{i-1}}{\Delta x^2}
\label{eq:WaveDisc1}
\end{equation}
\]
The above equation, \(\eqref{eq:WaveDisc1}\) is just a discretization of
Equation \(\eqref{eq:Wave1D}\) above, with
\(\frac{\partial^2 h}{\partial t^2}\) taken as acceleration. However, unlike the
simple spring system, where the position of the spring varied only in time,
our position varies in space as well, which means we have to consider that
every point in the array of heights depends on not only its previous positions
and velocities, but the heights and velocities and accelerations of neighboring
points. This was the oversimplification in our previous attempt, and the source
of the instability.
In order to understand this system better, let's replace the acceleration with a finite difference approximation, based on the height and previous heights. Starting with the backwards finite difference approximation of \(a_i\):
\[
\begin{equation}
a_{i_{t+\Delta t}} \approx \frac{h_{i_{t+\Delta t}} -
\ 2 h_{i_{t}} +
\ h_{i_{t-\Delta t}}}{\Delta t^2}
\label{eq:FiniteDiffA}
\end{equation}
\]
we can follow the usual convention in height-field solvers, eliminating
acceleration
by substituting Equation \(\eqref{eq:FiniteDiffA}\) into Equation
\(\eqref{eq:WaveDisc1}\) to get the following expression which contains
only height values (at different times and positions):
\[
\begin{equation}
\frac{h_{i_{t+\Delta t}} -
\ 2 h_{i_{t}} +
\ h_{i_{t-\Delta t}}}{\Delta t^2}
=
c^2 \frac{h_{{i+1}_{t+\Delta t}} -
\ 2 h_{i_{t+\Delta t}} +
\ h_{{i-1}_{t+\Delta t}} }
{\Delta x^2}
\label{eq:WaveHeightDisc1}
\end{equation}
\]
However, our RK4 solver and our other time-integration techniques were written
to employ a calculation of acceleration based on an estimate of current
position and velocity, so it's more helpful if we instead begin with the
following implicit relationship between acceleration, velocity, and height:
\[
\begin{equation}
h_{i_{t+\Delta t}} = h_{i_t} + \Delta t\ v_{i_{t+\Delta t}} +
\Delta t^2\ a_{i_{t+\Delta t}}
\label{eq:ImplicitHeightTimeStep}
\end{equation}
\]
Let's treat the first two terms of the right hand side of Equation
\(\eqref{eq:ImplicitHeightTimeStep}\) as an
estimate of height at time \(t+\Delta t\), denoted by \(h^\star\),
indicated as follows:
\[
\begin{equation}
h^\star_{i_{t+\Delta t}} = h_{i_t} + \Delta t\ v_{i_{t+\Delta t}}
\label{eq:HeightEstimate}
\end{equation}
\]
And then we can substitute Equation \(\eqref{eq:HeightEstimate}\) into
Equation \(\eqref{eq:ImplicitHeightTimeStep}\) to get:
\[
\begin{equation}
h_{i_{t+\Delta t}} = h^\star_{i_{t+\Delta t}} +
\Delta t^2\ a_{i_{t+\Delta t}}
\label{eq:EstimatedHeightTimeStep}
\end{equation}
\]
Now, we can subsitute Equation \(\eqref{eq:EstimatedHeightTimeStep}\) into
Equation \(\eqref{eq:WaveDisc1}\) to get the following equation, written
in terms of the current acceleration and an estimate of the current
height.
\[
\begin{equation}
a_{i_{t+\Delta t}} =
c^2 \frac{h^\star_{{i+1}_{t+\Delta t}} -
\ 2 h^\star_{i_{t+\Delta t}} +
\ h^\star_{{i-1}_{t+\Delta t}} }
{\Delta x^2} +
c^2 \Delta t^2
\frac{a_{{i+1}_{t+\Delta t}} -
\ 2 a_{i_{t+\Delta t}} +
\ a_{{i-1}_{t+\Delta t}} }
{\Delta x^2}
\label{eq:AccelDisc1Subscripted}
\end{equation}
\]
All of the terms in Equation \(\eqref{eq:AccelDisc1Subscripted}\) are at the
same point in time, so we can drop the time subscript. Assuming that the
estimates of height are given by some process (which we are free to tinker with),
the equation is a linear relationship between accelerations at different points
in space. Let's
simplify by defining intermediate constants:
\[
\begin{equation}
\kappa = \frac{c^2 \Delta t^2}{\Delta x^2},
\ \gamma = \frac{c^2}{\Delta x^2}
\label{eq:Constants}
\end{equation}
\]
Substituting in our constants, moving all of the acceleration terms to the
left-hand side, and gathering coefficients, we have the following equation
which expresses the relationship between the acceleration at the index \(i\)
and its spatially adjacent neighbors:
\[
\begin{equation}
(1 + 2 \kappa) a_i +
( -\kappa ) a_{i-1} +
( -\kappa ) a_{i+1} =
\gamma ( h^\star_{i+1} -
\ 2 h^\star_i +
\ h^\star_{i-1} )
\label{eq:AccelOneMatrixRow}
\end{equation}
\]
The above Equation \(\eqref{eq:AccelOneMatrixRow}\) is written relative to
any position in the array of heights, denoted by the subscript variable \(i\),
and this relationship exists at each point in the array.
We'll next write out the equations at each index explicitly, to get a system of
linear equations. We have to take a bit of care at the boundary points.
Our boundary condition
is that heights at the boundary are equal to the height at the adjacent,
non-boundary position, which, when transformed to a boundary condition on
acceleration, says that the acceleration at the boundary is zero. Therefore, we
have no equation for indices \(i = 0\) and \(i = N-1\), because we've
explicitly defined those points.
Furthermore, for indices \(i = 1\) and \(i = N-2\),
the relationship is slightly changed, which we'll incorporate into the
equations.
\[
\begin{eqnarray*}
\ &\ & (1+2\kappa)a_1& +& (-\kappa)a_2& =&
\gamma ( - h^\star_1 + h^\star_2 ) \\
(-\kappa)a_1& +& (1+2\kappa)a_2& +& (-\kappa)a_3& =&
\gamma ( h^\star_1 - 2 h^\star_2 + h^\star_3 ) \\
(-\kappa)a_2& +& (1+2\kappa)a_3& +& (-\kappa)a_4& =&
\gamma ( h^\star_2 - 2 h^\star_3 + h^\star_4 ) \\
... \\
(-\kappa)a_{n-5}& +& (1+2\kappa)a_{n-4}& +& (-\kappa)a_{n-3}& =&
\gamma ( h^\star_{n-5} - 2 h^\star_{n-4} + h^\star_{n-3} ) \\
(-\kappa)a_{n-4}& +& (1+2\kappa)a_{n-3}& +& (-\kappa)a_{n-2}& =&
\gamma ( h^\star_{n-4} - 2 h^\star_{n-3} + h^\star_{n-2} ) \\
(-\kappa)a_{n-3}& +& (1+2\kappa)a_{n-2}&\ &\ & =&
\gamma ( h^\star_{n-3} - h^\star_{n-2} ) \\
\label{eq:AccelLinearSystem}
\end{eqnarray*}
\]
This above system of linear equations can be expressed as a sparse matrix
equation:
\[
\begin{equation}
A x = b
\end{equation}
\]
Where the sparse, symmetric matrix \(A\) is defined as:
\[
\begin{equation}
A =
\begin{bmatrix}
1+2\kappa & -\kappa & \ & \ & \ & \ & \cdots \\
-\kappa & 1+2\kappa & -\kappa & \ & \ & \ & \cdots \\
\ & -\kappa & 1+2\kappa & -\kappa & \ & \ & \cdots \\
\ \vdots & \ & \ddots & \ddots & \ddots & \ & \cdots \\
\ \cdots & \ & \ & -\kappa & 1+2\kappa & -\kappa & \ \\
\ \cdots & \ & \ & \ & -\kappa & 1+2\kappa & -\kappa \\
\ \cdots & \ & \ & \ & \ & -\kappa & 1+2\kappa \\
\end{bmatrix}
\end{equation}
\]
and the column vectors \(x\) and \(b\) are defined as:
\[
\begin{equation}
x =
\begin{bmatrix}
a_1 \\
a_2 \\
a_3 \\
\vdots \\
a_{n-4} \\
a_{n-3} \\
a_{n-2} \\
\end{bmatrix}
,
b =
\begin{bmatrix}
\gamma ( - h^\star_1 + h^\star_2 ) \\
\gamma ( h^\star_1 - 2 h^\star_2 + h^\star_3 ) \\
\gamma ( h^\star_2 - 2 h^\star_3 + h^\star_4 ) \\
\vdots \\
\gamma ( h^\star_{n-5} - 2 h^\star_{n-4} + h^\star_{n-3} ) \\
\gamma ( h^\star_{n-4} - 2 h^\star_{n-3} + h^\star_{n-2} ) \\
\gamma ( h^\star_{n-3} - h^\star_{n-2} ) \\
\end{bmatrix}
\end{equation}
\]
We have thus transformed our acceleration computation problem into
an \(A x = b\) matrix problem, a problem which has been carefully studied,
and for which many extremely well-understood methods of solution exist.
We'll explore simple approaches to solving this system in the next post.
No comments:
Post a Comment
Note: Only a member of this blog may post a comment.