# Discretization for PDEs by Sn69ct4Y

VIEWS: 17 PAGES: 38

• pg 1
```									Discretization for PDEs

Chunfang Chen
Danny Thorne
Classification of PDEs
Different mathematical and physical
behaviors:
 Elliptic Type
 Parabolic Type
 Hyperbolic Type

System of coupled equations for several
variables:
 Time : first-derivative (second-derivative for
wave equation)
 Space: first- and second-derivatives
Classification of PDEs (cont.)
General form of second-order PDEs ( 2
variables)
PDE Model Problems
 Hyperbolic (Propagation)

• Wave equation    (Second-order linear )
PDE Model Problems (cont.)
 Parabolic (Time- or space-marching)

• Burger’s equation(Second-order nonlinear)

(Diffusion / dispersion)

• Fourier equation (Second-order linear )
PDE Model Problems (cont.)
 Elliptic (Diffusion, equilibrium problems)

• Laplace/Poisson (second-order linear)

• Helmholtz equation (second-order linear)
PDE Model Problems (cont.)
 System of Coupled PDEs

• Navier-Stokes Equations
Well-Posed Problem
Numerically well-posed
 Discretization equations
 Auxiliary conditions (discretized approximated)

• the computational solution exists (existence)
• the computational solution is unique (uniqueness)
• the computational solution depends continuously on the
approximate auxiliary data
• the algorithm should be well-posed (stable) also
Boundary and Initial
Conditions
Initial conditions: starting         R
point for propagation problems
R
Boundary conditions: specified
on domain boundaries to provide
the interior solution in          s
computational domain
n
Numerical Methods
 Complex geometry
 Complex equations (nonlinear, coupled)
 Complex initial / boundary conditions

 No analytic solutions
 Numerical methods needed !!
Numerical Methods
Objective: Speed, Accuracy at minimum
cost
 Numerical   Accuracy (error analysis)
 Numerical   Stability (stability analysis)
 Numerical   Efficiency (minimize cost)
 Validation  (model/prototype data, field data,
analytic solution, theory, asymptotic solution)
 Reliability and Flexibility (reduce preparation
and debugging time)
 Flow Visualization (graphics and animations)
computational solution
procedures
Governing                        System of      Equation    Approximate
Equations      Discretization    Algebraic      (Matrix)     Solution
ICS/BCS                          Equations       Solver

Continuous   Finite-Difference   Discrete    Tridiagonal    Ui (x,y,z,t)
Solutions                        Nodal
Values
Finite-Element                    SOR
T (x,y,z,t)
Spectral                    Gauss-Seidel
or
Boundary Element                  Krylov
Hybrid                      Multigrid      (,,, )
DAE
Discretization
 Time derivatives
almost exclusively by finite-difference methods
 Spatial derivatives
- Finite-difference: Taylor-series expansion
- Finite-element: low-order shape function and
interpolation function, continuous within each
element
- Finite-volume: integral form of PDE in each
control volume
- There are also other methods, e.g. collocation,
spectral method, spectral element, panel
method, boundary element method
Finite Difference
 Taylor series

 Truncation error

 How to reduce truncation errors?
• Reduce grid spacing, use smaller x = x-xo
• Increase order of accuracy, use larger n
Finite Difference Scheme

 Forward difference

 Backward difference

 Central difference
Example : Poisson Equation
(-1,1)    (1,1)

y

x

(-1,-1)   (1,-1)
Example (cont.)
Ui    1, j    2Ui , j  Ui     1, j       Ui ,   j  1    2Ui , j  Ui ,   j  1
                                             1
2                                             2
x                                            y
                                   
             (i,j+1)y             
     (i-1,j)x (i,j)(i+1,j)        
             (i,j-1)               
                                   
Rectangular Grid
After we discretize the Poisson equation on a
rectangular domain, we are left with a finite
number of gird points. The boundary values
of the equation are
0,4 1,4 2,4 3,4 4,4
the only known grid
0,3 1,3 2,3 3,3 4,3
points
0,2   1,2   2,2   3,2   4,2
0,1   1,1   2,1   3,1   4,1

0,0   1,0   2,0   3,0   4,0
What to solve?
Discretization produces a linear system of
equations.             4 1       0    1   0
The A matrix is a          1    4    1    0    1
pentadiagonal banded       0    1    4    1    0
matrix of a standard
1     0    1    4   1
form:
0    1     0    1    4
A solution method is to be performed for
solving
Matrix Storage
We could try and take advantage of the
banded nature of the system, but a
more general solution is the adoption of
a sparse matrix storage strategy.
Limitations of Finite
Differences
 Unfortunately, it is not easy to use
finite differences in complex
geometries.
 While it is possible to formulate
curvilinear finite difference methods,
the resulting equations are usually
pretty nasty.
Finite Element Method
 The finite element method, while more
complicated than finite difference methods,
easily extends to complex geometries.
 A simple (and short) description of the finite
element method is not easy to give.

Weak                Matrix
PDE
Form               System
Finite Element Method
(Variational Formulations)
 Find u in test space H such that a(u,v) = f(v) for all v
in H, where a is a bilinear form and f is a linear
functional.

     V(x,y) = j Vjj(x,y), j = 1,…,n
     I(V) = .5 j j AijViVj - j biVi, i,j = 1,…,n
     Aij =  a( j,  j), i = 1,…,n
     Bi =  f j, i = 1,…,n

 The coefficients Vj are computed and the function
V(x,y) is evaluated anyplace that a value is needed.
 The basis functions should have local support (i.e.,
have a limited area where they are nonzero).
Time Stepping Methods
 Standard methods are common:
– Forward Euler (explicit)
– Backward Euler (implicit)
– Crank-Nicolson (implicit)

T jn 1  T jn                                             θ = 0, Fully-Explicit
  LxxT    n 1
 (1   ) LxxT    n

t
j                         j
θ = 1, Fully-Implicit
T j 1  2T j  T j 1
LxxT j                                            θ = ½, Crank-Nicolson
x 2
Time Stepping Methods (cont.)
 Variable length time stepping
– Most common in Method of Lines (MOL) codes or
Differential Algebraic Equation (DAE) solvers
Solving the System
 The system may be solved using simple
iterative methods - Jacobi, Gauss-Seidel,
SOR, etc.

- No explicit storage of the matrix is required
- The methods are fairly robust and reliable

- Really slow (Gauss-Seidel)
- Really slow (Jacobi)
Solving the System
CG is a much more powerful way to solve the
problem.
– Easy to program (compared to other advanced
methods)
– Fast (theoretical convergence in N steps for an N
by N system)
– Explicit representation of the matrix is probably
necessary
– Applies only to SPD matrices
Multigrid Algorithm:
Components
 Residual
compute the error of the approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the
original grid
Residual Vector
 The equation we are to solve is defined as:

 So then the residual is defined to be:

 Where uq is a vector approximation for u
 As the u approximation becomes better,
the components of the residual vector(r) ,
move toward zero
Multigrid Algorithm:
Components
 Residual
compute the error of your approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the original
grid
Multigrid Algorithm:
Components
 Residual
compute the error of your approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the original
grid
The Restriction Operator
 Defined as ‘half-weighted’ restriction.
 Each new point in the courser grid, is
dependent upon it’s neighboring points from
the fine grid
Multigrid Algorithm:
Components
 Residual
compute the error of your approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the original
grid
The Prolongation Operator
 The grid change is exactly the opposite
of restriction

(0,2)    (1,2)    (2,2)    (0,4) (1,4) (2,4) (3,4) (4,4)
(0,3) (1,3) (2,3) (3,3) (4,3)

(0,1)    (1,1)    (2,1)    (0,2) (1,2) (2,2) (3,2) (4,2)
(0,1) (1,1) (2,1) (3,1) (4,1)
(0,0)    (1,0)    (2,0)    (0,0) (1,0) (2,0) (3,0) (4,0)
Prolongation vs. Restriction
 The most efficient multigrid algorithms use
prolongation and restriction operators that
are directly related to each other. In the one
dimensional case, the relation between
prolongation and restriction is as follows:
Full Multigrid Algorithm

1.Smooth initial U vector to receive a new approximation Uq
2. Form residual vector: Rq =b -A Uq
3. Restrict Rq to the next courser grid  Rq-1
4. Smooth Ae= Rq-1 starting with e=0 to obtain eq-1
5.Form a new residual vector using: Rq-1= Rq-1 -A eq-1
6. Restrict R2 (5x? where ?5) down to R1(3x? where ?3)
7. Solve exactly for Ae= R1 to obtain e1
8. Prolongate e1e2 & add to e2 you got from restriction
9. Smooth Ae= R2 using e2 to obtain a new e2
10. Prolongate eq-1 to eq and add to Uq
Reference

www.mgnet.org/
http://csep1.phy.ornl.gov/CSEP/PDE/PDE.html
www.ceprofs.tamu.edu/hchen/
www.cs.cmu.edu/~ph/859B/www/notes/multigrid.pdf
www.cs.ucsd.edu/users/carter/260
www.cs.uh.edu/~chapman/teachpubs/slides04-
methods.ppt
http://www.ccs.uky.edu/~douglas/Classes/cs521-
s01/index.html
Homework
Introduce the following processing
by read the book “A Tutorial on Elliptic PDE
Solvers and Their Parallelization”

Weak             Discrete
PDE
Form              Matrix

```
To top