UNIVERSITY OF WASHINGTON Department of Aeronautics and Astronautics

105 Pages·2003·3.46 MB·English

UNIVERSITY OF WASHINGTON
Department of Aeronautics and Astronautics

Numerical Techniques of Partial
Differential Equations

Prepared for Professor Nathan Kutz

by

Christopher Lum
Dec. 3, 2003
Table of Contents

Abstract...............................................................................................................................3
Introduction.........................................................................................................................3
Acknowledgments..............................................................................................................3
Nomenclature......................................................................................................................4
Theory.................................................................................................................................5
Finite Difference Discretization..................................................................................5
Elliptic Solve Using Fast Fourier Transform............................................................11
Elliptic Solve Using A\b...........................................................................................12
Elliptic Solve Using LU Decomposition..................................................................12
Elliptic Solve Using GMRES, BICGSTAB.............................................................13
Vorticity Coefficient.................................................................................................14
Implementation.................................................................................................................16
Control Flow.............................................................................................................16
Computation Difficulties..........................................................................................18
Results...............................................................................................................................23
Fast Fourier Transform.............................................................................................23
A\b.............................................................................................................................27
LU Decomposition....................................................................................................28
BICGSTAB...............................................................................................................29
GMRES.....................................................................................................................32
Spectral Implementation...........................................................................................35
Conclusion........................................................................................................................62

Appendix A: MatLab Functions......................................................................................63
Appendix B: MatLab Code..............................................................................................65
main.m......................................................................................................................65
user_preferences.m...................................................................................................75
create_matrices.m.....................................................................................................79
create_vorticity.m.....................................................................................................80
advection_fft.m.........................................................................................................83
analyze_wave_function.m........................................................................................85
ask_question.m..........................................................................................................85
display_parameters.m...............................................................................................87
view_movie.m...........................................................................................................88
view_subplots.m.......................................................................................................90
vorticity_coefficient.m..............................................................................................92
Appendix C: Derivation of Equations.............................................................................93
Equations of Motion.................................................................................................93
Finite Difference Discretization................................................................................96
Elliptic Solve Using Fast Fourier Transform..........................................................102
Appendix D: Code User’s Manual.................................................................................104
Appendix E: Problem Statement....................................................................................105

2
Abstract

In many applications of science and engineering, it becomes necessary to solve
partial differential equations of varying complexity. Even for simpler systems, many of
these partial differential equations do not have analytical solutions. In these situations it
becomes necessary to rely on numerical techniques to determine the evolution of the
system. This report investigates several methods to discretize and solve these partial
differential equations by writing them as a large system of ordinary differential equations.
Introduction

Fluid flow is one of the most complex applications to partial differential
equations. Equations that are derived from conservation of mass, momentum, and energy
are often very complex and require very advanced numerical techniques to solve. Even
simplified versions of these equations still require numerical techniques to solve.
However, the benefits of being able to simulate fluid flow are invaluable. This report
analyzes a simple evolution of vortices in a shallow fluid over time. Although this is a
fairly simple simulation, it serves as a basis for more complex computational fluid
dynamic algorithms.

Acknowledgments

The author would like to thank many of the people that contributed to the results
contained in this report. First and foremost, the guidance of Professor Nathan Kutz of the
University of Washington was invaluable. Many hours were spent in his class and office
getting the code to work. In addition, many conversations with classmates and
colleagues shed light on many confusing subjects.

3
Nomenclature

Table 1: Nomenclature used in code and in report
Report MatLab Description
Name Name
A A matrix representation of Laplacian operator. ∇2
B B matrix representation of partial derivative w.r.t x
C C matrix representation of partial derivative w.r.t y
C vorticity coefficient
ω
∇ ∇=∂/∂x+∂/∂y
M number of points in x direction used for vorticity coefficient
L
m m number of discretizations in x and y direction.
N number of points in y direction used for vorticity coefficient
L
n n n = m2.
n number of derivatives. ie. ∂n /∂xn
O of order(). ie. O(n2)
k wave number for Fourier transform
k KX wave number in x direction
x
k KY wave number if y direction
y
L Lmatrix lower triangular matrix decomposition of A
L L computational domain (20)
P P permutation matrix from LU decomposition.
ψ psi streamfunction
ψˆ(cid:4) psi_t streamfunction transformed in both x and y direction
ψ psi_col discretized streamfunction represented as a column vector
ψ ∂2ψ/∂x2. Similar notation for derivatives in different directions.
xx
t time required to solve elliptic equation
required
U Umatrix upper triangular matrix decomposition of A
υ v diffusion coefficient (0.001)
ω w vorticity
ω w_col discretized vorticity represented as a column vector
ωˆ(cid:4) w_t vorticity transformed in both x and y direction
ω wo initial vorticity
o
w.r.t with respect to
X Length of computational domain in x direction (equal to L here)
L
x x direction
y y direction
Y length of computational domain in y direction (equal to L here)
L

4
1
Theory

Presented here is a short summary of the general theory that governs this problem.
For a complete derivation of these equations and the procedure, please see Appendix C:
Derivation of Equations.

Finite Difference Discretization

The two governing equations of motion for this system are known as the elliptic
equation and the advection diffusion equation. These are given by Equation 1 and
Equation 2, respectively.

∇2ψ=ω
Equation 1

∂ω ∂ψ∂ω ∂ψ∂ω
=ν∇2ω+ −
∂t ∂y ∂x ∂x ∂y
Equation 2

As can be seen, these are partial differential equations in both the x and y
direction. Also, they are of varying order (the Laplacian is the sum of 2nd derivatives).
These are continuous operators and therefore, in order to numerically solve this problem,
these continuous operators must be discretized. Recall that the second order
discretization of a second derivative of a function is given by

d2f (t) f (t+∆t)−2f (t)+ f (t−∆t)
=
dt2 ∆t2
Equation 3

For convenience, we’ll choose ∆x=∆y =δ and for a given t, denote

ψ =ψ(m,n,t )
mn o

1 Reference [1]: Kutz, Nathan. PhD. University of Washington. Seattle, WA. Class notes for
AMATH581: Practical Scientific Computing. Autumn 2003.
5
With these simplifications, the elliptic equation (Equation 1) can be discretized to
the following form.

−4ψ +ψ +ψ +ψ +ψ
m,n m−1,n m+1,n m,n+1 m,n−1 =ω
δ2 m,n
Equation 4

From Equation 4, it can be seen that the vorticity at any point in the discretization
depends on the stream function at 4 neighboring points. Furthermore, with the periodic
boundary conditions

ψ =ψ
N+1,n 1,n
ψ =ψ
m,N+1 m,1

A similar procedure can be performed for discretizing the partial derivative of a
function in the x direction and the partial derivative of a function in the y direction. This
is shown below (recall that ∆x=∆y =δ).

∂ψ ψ −ψ
m,n = m+1,n m−1,n
∂x 2δ
Equation 5

∂ψ ψ −ψ
m,n = m,n+1 m,n−1
∂y 2δ
Equation 6

A simple stencil of this pattern for N = 4 is shown below in Figure 1

6
Figure 1: Pattern for vorticity and stream function at discrete points

2 matrices containing the x positions and the y positions can be computed using
the MESHGRID command

[X,Y] = meshgrid(x,y)

The initial vorticity over this grid can be generated using the X and Y matrices
and a given function.

ω = f (X,Y)
o

The important thing to notice is that the resulting matrix, ω, will have values of
o
vorticity at the different points stored in the following fashion.

ω ω ω ω 
11 21 31 41
 
ω ω ω ω
ω = 12 22 32 42
o ω ω ω ω 
13 23 33 43
 
 
ω ω ω ω
 
14 24 34 44
7
By comparing this with the stencil shown in Figure 1, one can realize that MatLab
stores values of vorticity that are reflected over the horizontal middle line. That is, the
value of vorticity in the upper left corner of the matrix ω is actually the value of
o
vorticity in the bottom left corner of the stencil. Similarly, the value of vorticity in the
upper right corner of the matrix ω is actually the value of vorticity in the bottom right
o
corner of the stencil. This is important because when the reshape command is used to
turn this matrix into a vector, the resulting column vector will be of the form.

wo_col = reshape(wo,16,1)

ω =(ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω )T
ocol 11 12 13 14 21 22 23 24 31 32 33 34 41 42 43 44

Therefore, when creating the discretization of the Laplacian operator (using
Equation 4) and the partial derivatives (using Equation 5 and Equation 6), the matrix
must be constructed such that it is compliant with this ordering of points. The Laplacian
can written in matrix form with this ordering is shown below as

−4 1 1 1 1 ψ  ω 
11 11
    
1 −4 1 1 1 ψ ω
  12  12
 1 −4 1 1 1 ψ  ω 
13 13
    
1 1 −4 1 1 1 ψ ω
  14  14
 1 1 −4 1 1 1 ψ  ω 
  21  21
 1 1 −4 1 1 ψ  ω 
22 22
    
1 1 −4 1 1 ψ ω
  23  23
1  1 1 1 −4 1 1 ψ  ω 
24 = 24
    
δ2 1 1 −4 1 1 1 ψ ω
  31  31
 1 1 −4 1 1 ψ  ω 
32 32
    
1 1 −4 1 1 ψ ω
  33  33
 1 1 1 −4 1 1 ψ  ω 
  34  34
 1 1 1 −4 1 1 ψ  ω 
41 41
    
1 1 1 −4 1 ψ ω
  42  42
 1 1 1 −4 1 ψ  ω 
43 43
    
 1 1 1 1 −4ψ  ω 
44 44

Equation 7

As can be seen, this is a sparse matrix with only 9 diagonals that result in
anything other than a zero entry. The large, sparse matrix can be referred to as the A
matrix, which performs the Laplacian operation on the vector representing the stream
8
function. This now becomes a large matrix equation of the form Aψ=ω. That is, given
a vorticity in column form, one can use a multitude of method to solve this for the stream
function vector.

The same procedure can be used to write the partial derivative in the x direction in
terms of a matrix equation. By comparing Equation 4 (Laplacian) with Equation 5
(partial w.r.t. x), one can see that the equation for the partial derivative is a simplified
version of the Laplacian operator. Therefore, the resulting matrix for the partial
derivative w.r.t x should be the same as the A matrix with less diagonal elements and
some with –1 instead of 1, and it should be divided by 2δ instead of δ2. The two matrix
operations which represent the partial w.r.t x and y, respectively, are shown below as
Equation 8 and Equation 9, respectively.

 1 −1 ψ  ψ 
11 11
    
1 −1 ψ ψ
  12  12
 1 −1 ψ  ψ 
13 13
    
1 −1 ψ ψ
  14  14
−1 1 ψ  ψ 
  21  21
 −1 1 ψ  ψ 
22 22
    
−1 1 ψ ψ
  23  23
1  −1 1 ψ  ∂ ψ 
24 = 24
    
2δ −1 1 ψ ∂x ψ
  31  31
 −1 1 ψ  ψ 
32 32
    
−1 1 ψ ψ
  33  33
 −1 1 ψ  ψ 
  34  34
 1 −1 ψ  ψ 
41 41
    
1 −1 ψ ψ
  42  42
 1 −1 ψ  ψ 
43 43
    
 1 −1 ψ  ψ 
44 44

Equation 8

9
 1 −1 ψ  ψ 
11 11
    
−1 1 ψ ψ
  12  12
 −1 1 ψ  ψ 
13 13
    
1 −1 ψ ψ
  14  14
 1 −1 ψ  ψ 
  21  21
 −1 1 ψ  ψ 
22 22
    
−1 1 ψ ψ
  23  23
1  1 −1 ψ  ∂ ψ 
24 = 24
    
2δ 1 −1 ψ ∂y ψ
  31  31
 −1 1 ψ  ψ 
32 32
    
−1 1 ψ ψ
  33  33
 1 −1 ψ  ψ 
  34  34
 1 −1ψ  ψ 
41 41
    
−1 1 ψ ψ
  42  42
 −1 1 ψ  ψ 
43 43
    
 1 −1 ψ  ψ 
44 44

Equation 9

The large, sparse matrices can be called B and C, which perform the operation of
partial derivative w.r.t x and partial derivative w.r.t y, respectively.

Using these matrices, the governing equation of motion for the evolution of the
fluid vorticity and stream function, Equation 1and Equation 2, can be rewritten as

Aψ=ω
Equation 10

∂ω
=CψBω−BψCω+υAω
∂t
Equation 11

A multitude of methods can be used to solve the elliptic problem (Equation 10)
for the stream function given a vorticity and a standard 2nd or 4th stepping scheme (ode23
or ode45) can be used to advance the advection-diffusion equation (Equation 11) into the
future. Once the new value of vorticity is found, the stream function can be updated and
the process can be repeated.
10