SimClass/BlogWaveMatrixMethods.html at master · blackencino/SimClass · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
<!DOCTYPE html>
<html>
<head>
  <title>Encino Sim Class : Matrix Methods and the Wave Equation</title>

  <style type="text/css" title="currentStyle">
    @import "css/SimClass.css";
  </style>

<script src='http://processingjs.org/js/processing.min.js' type='text/javascript'/></script>
<!Install the MathJax stuff so we can show LaTeX>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  TeX: { equationNumbers: { autoNumber: "AMS" } }
});
</script>
<script type="text/javascript"
  src="https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</head>
<body>
<h2>The Wave Equation as a Linear System</h2>
In our
<a href="http://encinographic.blogspot.com/2013/05/simulation-class-wave-equation.html">
previous class</a>,
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.
<p>
What follows in this post is a step-by-step progression, starting with the
original <b>PDE</b> 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.
<p>
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
<a href="https://en.wikipedia.org/wiki/Wave_equation">1D Wave Equation</a>
is written as:
<div class="LaTexEquation">
    \[
    \begin{equation}
    \frac{\partial^2 h}{\partial t^2} = c^2 \frac{\partial^2 h}{\partial x^2}
    \label{eq:Wave1D}
    \end{equation}
    \]
</div>
<p>
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\).
<div class="LaTexEquation">
    \[
    \begin{equation}
    a_i = c^2 \frac{h_{i+1} -\ 2 h_i + \ h_{i-1}}{\Delta x^2}
    \label{eq:WaveDisc1}
    \end{equation}
    \]
</div>
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.
<p>
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\):
<div class="LaTexEquation">
    \[
    \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}
    \]
</div>
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):

<div class="LaTexEquation">
    \[
    \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}
    \]
</div>

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 <i>implicit</i> relationship between acceleration, velocity, and height:

<div class="LaTexEquation">
    \[
    \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}
    \]
</div>

Let's treat the first two terms of the right hand side of Equation
\(\eqref{eq:ImplicitHeightTimeStep}\) as an
<i>estimate</i> of height at time \(t+\Delta t\), denoted by \(h^\star\),
indicated as follows:

<div class="LaTexEquation">
    \[
    \begin{equation}
    h^\star_{i_{t+\Delta t}} = h_{i_t} + \Delta t\ v_{i_{t+\Delta t}}
    \label{eq:HeightEstimate}
    \end{equation}
    \]
</div>

And then we can substitute Equation \(\eqref{eq:HeightEstimate}\) into
Equation \(\eqref{eq:ImplicitHeightTimeStep}\) to get:

<div class="LaTexEquation">
    \[
    \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}
    \]
</div>

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 <i>estimate</i> of the current
height.

<div class="LaTexEquation">
    \[
    \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}
    \]
</div>

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:

<div class="LaTexEquation">
    \[
    \begin{equation}
    \kappa = \frac{c^2 \Delta t^2}{\Delta x^2},
    \ \gamma = \frac{c^2}{\Delta x^2}
    \label{eq:Constants}
    \end{equation}
    \]
</div>

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:

<div class="LaTexEquation">
    \[
    \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}
    \]
</div>

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.

<div class="LaTexEquation">
    \[
    \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*}
    \]
</div>

This above system of linear equations can be expressed as a sparse matrix
equation:
<div class="LaTexEquation">
    \[
    \begin{equation}
    A x = b
    \end{equation}
    \]
</div>
Where the sparse, symmetric matrix \(A\) is defined as:
<div class="LaTexEquation">
\[
\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}
 \]
</div>
and the column vectors \(x\) and \(b\) are defined as:
<div class="LaTexEquation">
\[
\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}
 \]
</div>
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.

<p>
<h2>Download</h2>
All of the files associated with this class are available via a public github repository, found here:
<a href="https://github.com/blackencino/SimClass">https://github.com/blackencino/SimClass</a>

</body>