7/20/2017

Instructor
Dr. Raymond Rumpf
(915) 747‐6958
rcrumpf@utep.edu

EE 5337
Computational Electromagnetics (CEM)

Lecture #10

Finite‐Difference Method
 These notes may contain copyrighted material obtained under fair use rules.  Distribution of these materials is strictly prohibited  
Lecture 10 Slide 1

Outline
• Finite‐Difference Approximations
• Finite‐Difference Method
• Numerical Boundary Conditions
• Matrix Operators

Lecture 10 Slide 2

1
 7/20/2017

Finite Difference 
Approximations

Lecture 10 Slide 3

The Basic Finite‐Difference Approximation

df1.5 f 2  f1
 second‐order accurate

x
first‐order derivative
dx
df
dx
f2
f1
x

This is the only finite‐difference approximation we will use in this course!
Lecture 10 Slide 4

2
 7/20/2017

Types of Finite‐Difference Approximations

df 2 f 2  f1
Backward 

Finite‐Difference
dx x

df1.5 f 2  f1
Central 

Finite‐Difference
dx x

df1 f 2  f1
Forward 

Finite‐Difference
dx x

Lecture 10 Slide 5

The Generalized Finite‐Difference

dn f
  ai fi L  f    ai fi
dx n i i
In fact, any linear operation on the 
function can be approximated as a 
linear sum of known points of that 
function. 

The derivative of any order of a 
function at any position can be 
approximated as a linear sum of 
known points of that function. 

Lecture 10 Slide 6

3
 7/20/2017

Finite‐Difference “Atoms”
Any finite‐difference approximation can be summarized graphically as 
an “atom.”
1 1
2
  2 d f f
f  xi   i 1 i 1
f i 1 fi f i 1 dx 2

1 2 1

2
 2  2
d2 f  2 fi  fi 1
f  xi   i 1
f i 1 fi f i 1 dx 2
2

1 2
f i , j 1
1  1
2 4 fi, j 2 f i 1, j  2 f i , j  fi 1, j f i , j 1  2 f i , j  fi , j 1

2
 f i 1, j  2 f  xi   
f i 1, j  2
2

1 2
f i , j 1
Lecture 10 Slide 7

Finite‐Difference
Method

Lecture 10 Slide 8

4
 7/20/2017

Overview of Our FDM
1. Identify and write the governing equation(s).
2 
a  x f  x   b  x f  x  c  x f  x  g  x
x 2 x
2. Write the matrix form of this equation going term‐by‐term.
2 
a  x f  x   b  x f  x  c  x f  x  g  x
x 2
x

 A  Dx2   f     B  Dx  f   C  f    g 
3. Put matrix equation in the standard form of [L][f] = [g].
 L  f    g  where  L    A  Dx2     B  Dx   C 
L = A*D2X + y*B*DX + C;
4. Solve [L][f] = [g].
 f    L  g 
1

Lecture 10 Slide 9

What’s the Catch?

 L    A  Dx2     B  Dx   C 

How do we construct these matrices?  What do they mean?
 A
 Dx 
  B  Dx2 
C 

Lecture 10 Slide 10

5
 7/20/2017

Functions Vs. Operations (1 of 2)

2 
a  x f  x   b  x f  x  c  x f  x  g  x
x 2 x

Operations Functions
Everything else in a differential  The only time functions 
equation is something that  appear in a differential 
operates on a function. equation is as the unknown 
or as the excitation.
a  x  , b  x  , c  x   point-by-point
multiplication on f  x  f  x   unknown
  2
,  calculates derivatives of f  x  g  x   excitation
x x 2
  scales entire f  x 

Lecture 10 Slide 11

Functions Vs. Operations (2 of 2)

 A  Dx2   f     B  Dx  f   C  f    g 

Operations Functions
Operations are always stored  Functions are stored as 
in square matrices.  Any  column vectors.
linear operation can be put 
into matrix form.  f1   g1 
f  g 
 l11 l12 l1M   f    2   g    2 
l l22  l2 M     
  
L   21

 
 fM   gM 
 
lM 1 lM 2 lMM 

Lecture 10 Slide 12

6
 7/20/2017

Interpretation of Matrices
a11 x  a12 y  a13 z  b1
 a11 a12 a13   x   b1 
a21 x  a22 y  a23 z  b2 a
 21 a22 a23   y   b2 
a31 x  a32 y  a33 z  b3  a31 a32 a33   z  b3 

EQUATION FOR… RELATION TO…

Equation for x  a11 a12 a13 

a a22 a23   a11 a12 a13 
Equation for y  21 a
Equation for z  a31 a32 a33   21 a22 a23 
 a31 a32 a33 
From a purely mathematical perspective, this interpretation does not make sense.  This 
interpretation will be highly useful and insightful because of how we derive the equations.
Lecture 10 Slide 13

Representing Functions on a Grid

Example  A grid is 
physical  constructed by  Function is  Representation 
(continuous)  dividing space  known only at  of what is 
2D function into discrete  discrete points actually stored in 
cells memory

Lecture 10 Slide 14

7
 7/20/2017

Grid Cells

A Single Grid Cell
x
Whole Grid
y

A function value is 
assigned to a specific 
point within the grid 
cell.

x, y  grid resolution parameters

Lecture 10 Slide 15

Functions are Put Into Column Vectors

1‐D Systems  f1 
f 
 2
2‐D Systems  f3 
 f1   
 f4 
f   f5 
f1
f 2 f3 f 4 f5  2
f   f3  f1 f5 f9 f13
 
 f6 
  f 
 f4  f2 f6 f10 f14  7
 f5  f 
f  8
f3 f7 f11 f15 f
 9
 f10 
 
f4 f8 f12 f16  f11 
 f12 
 
 f13 
f 
 14 
 f15 
f 
 16 

Lecture 10 Slide 16

8
 7/20/2017

Putting Functions into Column Vectors

 f1 
f1 f 
f2  2
 f3 
f3  
f4
 f4 
 f5 
f1 f5 f9 f13 f5  
 f6 
f6 f 
f2 f6 f10 f14  7
f7 f 
f  8
f3 f7 f11 f15 f8 f
 9
f9  f10 
f4 f8 f12 f16  
f10  f11 
f11  f12 
 
f12  f13 
f 
f13  14 
 f15 
f14 f 
MATLAB ‘reshape’ command f15  16 
f = F(:); f16
F = reshape(f,Nx,Ny);

Lecture 10 Slide 17

Locating Nodes in Column Vectors
1D Grids
Node located at nx  m = nx

2D Grids
Node located at nx,ny  m = (ny – 1)*Nx + nx

3D Grids
Node located at nx,ny,nz
 m = (nz – 1)*Nx*Ny + (ny – 1)*Nx + nx

Lecture 10 Slide 18

9
 7/20/2017

The Finite‐Difference Method
Differential Equation to Solve
df  x 
 a  x f  x  g  x
dx
Tedious, but 
Almost effortless.
Conventional FDM

not difficult.

Improved FDM
f  k  1  f  k  1
 a k  f k   g k   Dx  f    a  f    g 
2x

Very difficult  Very easy 
and tedious. and clean

 L  f    g 
 L  f    g 
 L    Dx    a 
Easy. Easy.

 f    L  g 
1

Most time consuming step.

Lecture 10 Slide 19

Conventional FDM (1 of 3)
Step 1 – We start with a differential equation that we wish to solve.
d2 f df
2
 a  bf  c
dx dx
Step 2 – We approximate the derivatives with finite differences.
f  k  1  2 f  k   f  k  1 f  k  1  f  k  1
 a k   b k  f k   c k 
 2
2
IMPORTANT RULE:  Every term in a finite‐difference equation must exist at the 
same point.

Step 3 – The equation is expanded and we collect common terms
1 2 1 a k  a k 
f  k  1  2 f  k   2 f  k  1  f  k  1  f  k  1  b  k  f  k   c  k 
2
   2 2

 1 a k    2   1 a k  
 2  f  k  1  b  k   2  f  k    2   f  k  1  c  k 
  2       2 
Lecture 10 Slide 20

10
 7/20/2017

Conventional FDM (2 of 3)
Step 4 – The final equation is used to populate a matrix equation.
 1 a k    2  1 a k  
 2  f  k  1  b  k   2  f  k    2   f  k  1  c  k 
 2      2 

  f 1   c 1 
  f  2
  
Filling in the      c  2 
matrix like    f  3   c  3 
this can be a      
very difficult    f  4   c  4 
and tedious    f  5   c  5 
task for more     
complicated      
differential        
equations or      
for systems of      
differential   f  N  1  c  N  1 
 
equations.  
   f  N    c  N  

Lecture 10 Slide 21

Conventional FDM (3 of 3)
Step 5 – The matrix equation is solved for the unknown function f(x)

 L  f   c 
Lf  c

 f    L c
1

f  L1c

Lecture 10 Slide 22

11
 7/20/2017

Improved FDM
We want a very easy way to immediately write differential equations 
in matrix form.  Starting with the same differential equation…

2 f f
 a  bf  c
x 2
x
We will develop a procedure by which this will be directly written in 
matrix form without having to explicitly handle any finite‐differences.

2  Dx  , A, Dx  , and B

2 1
f a f  bf  c
x x
2 are square matrices that perform linear 
operations on the vector f.

Dx2f  ADx1f  Bf = c  Dx2  ADx1  B  f = c


  
L  L 
f  L1c
Lecture 10 Slide 23

Locating Nodes in Matrices
Relation to node (nx,ny,nz)
Column n

1D    m or n = nx
2D    m or n = (ny–1)*Nx + nx
3D    m or n = (nz–1)*Nx*Ny + (ny–1)*Nx + nx

Row m
Equation for node (nx,ny,nz)

Lecture 10 Slide 24

12
 7/20/2017

Numerical
Boundary Conditions

Lecture 10 Slide 25

The Problem at the Boundaries
Suppose it is desired to solve the following differential equation on a 
33 grid.
 2 f  x, y   2 f  x, y  fi 1, j  2 f i , j  f i 1, j f i , j 1  2 fi , j  f i , j 1
  g  x, y     gi , j
x 2
y 2
x 2
y 2

f1,0 f 2,0 f 3,0 f 2,1  2 f1,1  f 0,1 f1,2  2 f1,1  f1,0

  g1,1
x 2 y 2
f 0,1 f 4,1 f 3,1  2 f 2,1  f1,1 f 2,2  2 f 2,1  f 2,0
  g 2,1
x 2 y 2
f 4,1  2 f 3,1  f 2,1 f 3,2  2 f 3,1  f 3,0
  g 3,1
f 0,2 f 4,2 x 2 y 2
f 2,2  2 f1,2  f 0,2 f1,3  2 f1,2  f1,1
  g1,2
x 2 y 2
f 3,2  2 f 2,2  f1,2 f 2,3  2 f 2,2  f 2,1
f 0,3 f 4,3   g 2,2
x 2 y 2
f 4,2  2 f 3,2  f 2,2 f 3,3  2 f 3,2  f 3,1
  g 3,2
x 2 y 2
f1,4 f 2,4 f 3,4
f 2,3  2 f1,3  f 0,3 f1,4  2 f1,3  f1,2
  g1,3
x 2 y 2
The terms in red exist outside of the grid. f 3,3  2 f 2,3  f1,3

f 2,4  2 f 2,3  f 2,2
 g 2,3
x 2 y 2
How this is handled is called a boundary f 4,3  2 f 3,3  f 2,3 f 3,4  2 f 3,3  f 3,2
  g 3,3
condition. x 2 y 2

Lecture 10 Slide 26

13
 7/20/2017

Dirichlet Boundary Conditions (1 of 2)
The simplest boundary condition is to assume that all values of f(x,y)
outside of the grid are zero.

f1,0 f 2,0 f 3,0 f 2,1  2 f1,1  0 f1,2  2 f1,1  0

  g1,1
x 2 y 2
f 0,1 f 4,1 f 3,1  2 f 2,1  f1,1 f 2,2  2 f 2,1  0
  g 2,1
x 2 y 2
0  2 f 3,1  f 2,1 f 3,2  2 f 3,1  0
  g 3,1
f 0,2 f 4,2 x 2 y 2
f 2,2  2 f1,2  0 f1,3  2 f1,2  f1,1
  g1,2
x 2 y 2
f 3,2  2 f 2,2  f1,2 f 2,3  2 f 2,2  f 2,1
f 0,3 f 4,3   g 2,2
x 2 y 2
0  2 f 3,2  f 2,2 f 3,3  2 f 3,2  f 3,1
  g 3,2
x 2 y 2
f1,4 f 2,4 f 3,4
f 2,3  2 f1,3  0 0  2 f1,3  f1,2
  g1,3
x 2 y 2
f 3,3  2 f 2,3  f1,3 0  2 f 2,3  f 2,2
  g 2,3
x 2 y 2
0  2 f 3,3  f 2,3 0  2 f 3,3  f 3,2
  g 3,3
x 2 y 2

Lecture 10 Slide 27

Dirichlet Boundary Conditions (2 of 2)
Dirichlet boundary conditions assume function values from outside 
of the grid are zero.

f  x1 , yi  f  x2 , yi   0

x 2 x

f xN x , yi   0 f x N x 1 , yi   2 f  x1 , yi  f  x2 , yi   2 f  x1 , yi   0

x 2 x x 2  2x
f  xi , y1  f  xi , y2   0 
 2 f xN x , yi   0 2 f x  
, yi  f xN x 1 , yi 
 Nx
y 2 y
x 2
 2
x


f xi , y N y  
0  f xi , yN y 1   f  xi , y1  f  xi , y2   2 f  xi , y1   0
2


y 2 y y 2  2y


 2 f xi , y N y   0 2 f x , y  f x , y  i Ny i N y 1

y 2
 2
y

Lecture 10 Slide 28

14
 7/20/2017

Periodic Boundary Conditions (1 of 2)
If the function f(x,y) is periodic, then the values from outside of the 
grid can be mapped to a value from inside the grid at the other side.
f 0, j  f 3, j , f 4, j  f1, j , fi ,0  f i ,3 , and f i ,4  fi ,1

f1,0 f 2,0 f 3,0 f 2,1  2 f1,1  f 3,1 f1,2  2 f1,1  f1,3

  g1,1
x 2 y 2
f 0,1 f 4,1 f 3,1  2 f 2,1  f1,1 f 2,2  2 f 2,1  f 2,3
  g 2,1
x 2 y 2
f1,1  2 f 3,1  f 2,1 f 3,2  2 f 3,1  f 3,3
  g 3,1
f 0,2 f 4,2 x 2 y 2
f 2,2  2 f1,2  f 3,2 f1,3  2 f1,2  f1,1
  g1,2
x 2 y 2
f 3,2  2 f 2,2  f1,2 f 2,3  2 f 2,2  f 2,1
f 0,3 f 4,3   g 2,2
x 2 y 2
f1,2  2 f 3,2  f 2,2 f 3,3  2 f 3,2  f 3,1
  g 3,2
x 2 y 2
f1,4 f 2,4 f 3,4
f 2,3  2 f1,3  f 3,3 f1,1  2 f1,3  f1,2
  g1,3
x 2 y 2
f 3,3  2 f 2,3  f1,3 f 2,1  2 f 2,3  f 2,2
  g 2,3
x 2 y 2
f1,3  2 f 3,3  f 2,3 f 3,1  2 f 3,3  f 3,2
  g 3,3
x 2 y 2

Lecture 10 Slide 29

Periodic Boundary Conditions (2 of 2)
Periodic boundary conditions assume function values from outside 
of the grid can be taken from the opposite side of the grid.

f  x1 , yi  f  x2 , yi   f xN x , yi

 
x 2 x

f xN x , yi   f x , y  f x
1 i N x 1 , yi   2 f  x1 , yi  f  x2 , yi   2 f  x1 , yi   f xN x , yi  
x 2 x 
x 2  2x
f  xi , y1  f  xi , y2   f xi , yN y

  
 2 f xN x , yi   f x , y 2 f x
1 i Nx  
, yi  f xN x 1 , yi 
y 2 y x 2
 2
x


f xi , y N y  
f  xi , y1   f xi , yN y 1   f  xi , y1  f  xi , y2   2 f  xi , y1   f xi , y N y
2
 
y 2 y 
y 2  2y


 2 f xi , y N y   f x , y 2 f x , y  f x , y 
i 1 i Ny i N y 1

y 2
 2
y

Lecture 10 Slide 30

15
 7/20/2017

Neumann Boundary Conditions
We use Neumann boundary conditions for non‐periodic functions 
or functions that are not zero at the boundary.  Here the function 
continues linearly outside of the grid.

Spatially variant grating  Here is a 1D function with 
that is not periodic and not  Neumann boundaries
zero at the boundaries.

Lecture 10 Slide 31

Neumann BC’s for 1D Function
The finite‐difference approximation for a 1D function is
df f f
 i 1 i 1
dx i 2 x

At i=1, we have a problem…

df f f This term  df f f
 2 0 doesn’t exist!  2 1
dx i 1 2 x dx i 1 x

At i=Nx, we have another problem…
This term 
df f N 1  f N x 1 doesn’t exist! df f N  f N x 1
 x  x
dx i  N x 2 x dx i  N x x
Lecture 10 Slide 32

16
 7/20/2017

What About the 2nd‐Order Derivatives for 
the Neumann Boundary Condition?
In order for the function to continue in a 
straight line, the second‐order derivate 
should be set to zero at the boundary.

IMPORTANT:  This is NOT Dirichlet boundary 
conditions.  A Dirichlet BC sets the function 
itself to zero outside of the grid, not the 
derivative.  Here, the 2nd‐order derivative is 
set to zero.

d2 f f f f d2 f
 2 12 0  0
dx 2 i 1
x dx 2 i 1

d2 f f N x 1  f N x  f N x 1 d2 f
  0
dx 2 i Nx
 2x dx 2 i Nx

Lecture 10 Slide 33

Neumann Boundary Conditions for 2D 
Functions
Neumann boundary conditions are used when a function should be 
continuous at the boundary.  That is, the first‐order derivative is 
continuous and the second‐order derivative is zero.

f  x1 , yi  f  x2 , yi   f  x1 , yi   2 f  x1 , yi 
 0
x x x 2

f xN x , yi   f x Nx  
, yi  f xN x 1 , yi  
 2 f xN x , yi  0
x x x 2

f  xi , y1  f  xi , y2   f  xi , y1   f  xi , y1 
2

 0
y y y 2


f xi , y N y   f x , y  f x , y 
i Ny i N y 1 
 2 f xi , yN y  0
y y y 2

Lecture 10 Slide 34

17
 7/20/2017

Matrix Operators

Lecture 10 Slide 35

Origin of Matrix Operators
We approximate the governing equation with 
We start with a governing  finite‐differences and then write the finite‐
equation. difference equation at each point the grid.
f 2,1  2 f1,1  0
 g1,1
x 2
f 3,1  2 f 2,1  f1,1
d 2 f  x, y  x 2
 g 2,1 We collect the large set of equations 
2
 g  x, y  0  2 f3,1  f 2,1 into a single matrix equation.
dx  g3,1
x 2
f 2,2  2 f1,2  0
 g1,2  2 1 0 0 0 0 0 0 0   f1,1   g1,1 
x 2    f  g 
f 3,2  2 f 2,2  f1,2
 g 2,2  1 2 1 0 0 0 0 0 0   2,1   2,1 
x 2  0 1 2 0 0 0 0 0 0   f3,1   g 3,1 
0  2 f 3,2  f 2,2      
 g3,2  0 0 0 2 1 0 0 0 0   f1,2   g1,2 
x 2 1     f 2,2    g 2,2 
f 2,3  2 f1,3  0
0 0 0 1  2 1 0 0 0
 g1,3 x 2      
x 2  0 0 0 0 1 2 0 0 0   f3,2   g 3,2 
f 3,3  2 f 2,3  f1,3  0 0 0 0 0 0 2 1 0  f    g 
 g 2,3    1,3   1,3 
x 2  0 0 0 0 0 0 1 2 1   2,3   2,3 
f g
0  2 f 3,3  f 2,3  0 0 0 0 0 0 0 1 2   f   g 
 g3,3    3,3   3,3 
x 2

We construct a grid to  D2x
store the functions. This matrix calculates the derivative of f(x,y)
and puts the answer in g(x,y).  This is a matrix 
Lecture 10
operator. Slide 36

18
 7/20/2017

Other Matrix Operators
A square matrix can always be constructed to perform any linear 
operation on a function that is stored in a column vector.

   f1 
  f 
d   2 
f  x  D xf   ?   f3 
dx   
   
   f N 

Dx

FFT  f  x    Ff g  x  f  x   Gf

  f  x  dx    f
 h  x   f  x   Hf

Lecture 10 Slide 37

Point‐by‐Point Multiplication (1 of 2)

b3 b4 b5 b6 Since we are storing our “functions” in vector form, 
b2 how do we perform a point‐by‐point multiplication 
b1 using a square matrix?

f1
f2 f3 f4
f5 f6
b  x f  x  Bf
x1 x2 x3 x4 x5 x6

b1 0 0 0 0 0   f1   b1 f1 
 0 b 0 0 0 0   f  b f 
 2  2  2 2

?
 0 0 b3 0 0 0   f 3   b3 f3 
   
 0 0 0 b4 0 0   f 4  b4 f 4 
 0 0 0 0 b5 0   f 5   b5 f5 
   



 f 6  b6 f 6 
 0 0 0 0 0 b6  
 
B f Bf
Lecture 10 Slide 38

19
 7/20/2017

Point‐by‐Point Multiplication (2 of 2)

b3 b4 b5 b6 Since we are storing our “functions” in vector form, 
b2 how do we perform a point‐by‐point multiplication 
b1 using a square matrix?

f1
f2 f3 f4
f5 f6
b  x f  x  Bf
x1 x2 x3 x4 x5 x6

b1 0 0 0 0 0   f1   b1 f1 
 0 b 0 0 0 0   f  b f 
 2  2  2 2
 0 0 b3 0 0 0   f 3   b3 f3 
     
 0 0 0 b4 0 0   f 4  b4 f 4 
 0 0 0 0 b5 0   f 5   b5 f5 
    

0 0 b6   f 6  b6 f 6 
 0 0 0    
B f Bf
Lecture 10 Slide 39

First‐Order Partial Derivative (1 of 2)
How do we construct a square matrix so that when it premultiplies a vector, we get 
a vector containing the first‐order partial derivative?


f  x  Dx f
1
f1
f6
f2 f3 f4
f5
x
x1 x2 x3 x4 x5 x6

 0 1 0 0 0 0   f1   f 2  f 0  2 x 
 1 0 1 0 0 0   f   f  f 2 
   2   3 1  x 

1 

0  1
?
0 1 0 0   f 3
   
   4 2 x
f  f 2
2  x  0 0 1 0 1 0   f 4    f 5  f 3  2  x 
 0 0 0 1 0 1   f 5    f 6  f 4  2  x 
   




 0 0  0 0 1 0   f 6   f 7  f 5  2 x 
   
Dx1 f Dx1f
Lecture 10 Slide 40

20
 7/20/2017

First‐Order Partial Derivative (2 of 2)
How do we construct a square matrix so that when it premultiplies a vector, we get 
a vector containing the first‐order partial derivative?


f  x  Dx f
1
f1
f6
f2 f3 f4
f5
x
x1 x2 x3 x4 x5 x6

 0 1 0 0 0 0   f1   f 2  f 0  2 x 
 1 0 1 0 0 0   f   f  f 2 
   2   3 1  x 

1  0 1 0 1 0 0   f 3    f 4  f 2  2  x 
     
2  x  0 0 1 0 1 0   f 4    f 5  f 3  2  x 
 0 0 0 1 0 1   f 5    f 6  f 4  2  x 
    
 0 0 0 0 1 0   f 6   f 7  f 5  2 x 
   

Dx1 f Dx1f
Lecture 10 Slide 41

Second‐Order Partial Derivative (1 of 2)
How do we construct a square matrix so that when it premultiplies a vector, we get 
a vector containing the second‐order partial derivative?

2
f1
f6 f  x  Dx2f
f2 f3 f4
f5
x 2
x1 x2 x3 x4 x5 x6

 2 1 0 0 0 0   f1    f 2  2 f1  f 0   x 
2

 1 2 1 0 0 0   f   f  2 f  f  2 
   2   3 2 1 x

1 

0 1 2
? 1 0 0  
 
f 3




 4
f  2 f 3  f 2

 2x  0 0 1 2 1 0   f 4   f 5  2 f 4  f3   2x 
 0 0 0 1 2 1   f 5   f 6  2 f 5  f 4   2x 
   
 2
x



 0 0 0 0 1 2   f 6   f 7  2 f 6  f 5   2x 
    
Dx2 f Dx2f
Lecture 10 Slide 42

21
 7/20/2017

Second‐Order Partial Derivative (2 of 2)
How do we construct a square matrix so that when it premultiplies a vector, we get 
a vector containing the second‐order partial derivative?

2
f1
f6 f  x  Dx2f
f2 f3 f4
f5
x 2
x1 x2 x3 x4 x5 x6

 2 1 0 0 0 0   f1    f 2  2 f1  f 0   x 
2

 1 2 1 0 0 0   f   f  2 f  f  2 
   2   3 2 1 x

1  0 1 2 1 0 0   f 3   f 4  2 f3  f 2   x 
2

      
 2x  0 0 1 2 1 0   f 4   f 5  2 f 4  f3   2x 
 0 0 0 1 2 1   f 5   f 6  2 f 5  f 4   2x 
    
 0 0 0 0 1 2   f 6   f 7  2 f 6  f 5   2x 
   

Dx2 f Dx2f
Lecture 10 Slide 43

What About 2D Grids (1 of 2)?
Two‐dimensional grids are a little more difficult x
f1 f2 f3

2
f  x, y   Dx f
2 f4 f5 f6

x 2 f7 f8 f9

y
 2 1   f1    f 2  2 f1  f ???   x 
2

 1 2 1   f   f  2 f  f 2 
  2   3 2 1 x 

 1 2 0   f 3
   f???  2 f3  f 2   2x 
    
 0 2 1   f 4   f 5  2 f 4  f ???   2x 
1    f5     f 6  2 f 5  f 4   2x 
1 2 1
 2x     
1 2 0   f 6   f???  2 f 6  f5   x 
2

 0 2 1   f   f  2 f  f   2 
  7  8 7 ??? x

 1 2 1   f8    f 9  2 f8  f 7   2x 
 

 1 2   f 9   f ???  2 f 9  f8   2x 
    
Dx2 f Dx2f
Lecture 10 Slide 44

22
 7/20/2017

What About 2D Grids (2 of 2)?
Two‐dimensional grids are a little more difficult x
f1 f2 f3

2
f  x, y   Dy2f
f4 f5 f6

y 2 f7 f8 f9

y
 2 0 0 1   f1    f 4  2 f1  f ???   y 
2

 0 2 0 0 1  f   f 2f  f 2
   2   5 2 ???   y 

 0 0 2 0 0 1   f 3    f 6  2 f3  f???   y 2

    
 1 0 0 2 0 0 1   f 4    f 7  2 f 4  f1   2y 
1    f 5     f8  2 f 5  f 2   2y 
1 0 0 2 0 0 1
 2y     
 1 0 0 2 0 0 1   f 6    f9  2 f 6  f3   2y 
 1 0 0 2 0 0   f 7   f ???  2 f 7  f 4   2y 
    
 1 0 0 2 0   f8    f???  2 f8  f 5   2y 
 

 1 0 0 2   f9   f ???  2 f9  f 6   2y 
    
Dy2 f Dy2f
Lecture 10 Slide 45

Derivative Operators with Dirichlet Boundary 
Conditions
 2 1 0 0 0 0 0 0 0  Both of these matrices only have 
 1 2 1 0 0 0 0 0 0  numbers along three of their diagonals.
 
 0 1 2 0 0 0 0 0 0  This is called tridiagonal.
 
0 0 0 2 1 0 0 0 0 
1 
This suggests a fast way to construct 
Dx2  2  0 0 0 1 2 1 0 0 0  these matrices.
x  
 0 0 0 0 1  2 0 0 0  Dx(2)  has some corrections shown in 
 0 0 0 0 0 0 2 1 0  blue along two of its diagonals.
 
 0 0 0 0 0 0 1 2 1 
 0 0 0 0 0 0 0 1 2   2 0 0 1 0 0 0 0 0 
 
 0 2 0 0 1 0 0 0 0 
 
These matrices contain mostly zeros.  0 0 2 0 0 1 0 0 0 
 
1  1 0 0 2 0 0 1 0 0 
 
These are called a sparse matrices. D y  2  0 1 0 0 2 0 0 1 0 
2

y  
See MATLAB sparse() command.  0 0 1 0 0 2 0 0 1 
 0 0 0 1 0 0 2 0 0 
 
Also see MATLAB spdiags() command.  0 0 0 0 1 0 0 2 0 
 0 0 0 0 0 1 0 0 2 
 
Lecture 10 Slide 46

23
 7/20/2017

Why Do We Need Separate Matrix Operators for 
First‐ and Second‐Order Derivatives?
We know that,
2   2  
 
x 2
x x y 2
y y
Can we just calculate D(2) from D(1)?
? ?
Dx   Dx  Dx  Dy   Dy  Dy 
2 1 1 2 1 1

Yes, but this does not make efficient use of the grid.  For a 5‐point, 
1D grid, we have
 1 0 1 0 0   2 1 0 0 0 
 0 2 0 1 0   1 2 1 0 0 
1   1  
Dx  Dx    1 0 2 0 1  Dx   2  0 1 2 1 0 
1 1 2

 2 x  x 
2
  
 0 1 0 2 0   0 0 1 2 1 
 0 0 1 0 1  0 0 0 1 2 
This is not as accurate because it calculates the  This matrix operator makes optimal use of 
derivative with poorer grid resolution than is available. the available grid resolution.
Lecture 10 Slide 47

2D Derivative Operators for 1D Grids
When Nx=1 and Ny>1
0 0
Dx  Z zero matrix  0 
Dx   
Dy is standard for 1D grid   
 
0 0

When Nx>1 and Ny=1
Dx is standard for 1D grid
0 0
Dy  Z zero matrix  0 
Dy   
  
 
0 0

Lecture 10 Slide 48

24
 7/20/2017

USE SPARSE MATRICES!!!!!!!

WARNING !!
The derivative operators will be EXTREMELY large matrices.
For a small grid that is just 100200 points:

Total Number of Points: 20,000

Size of Derivate Operators: 20,000  20,000
Total Elements in Matrices: 400,000,000
Memory to Store One Full Matrix: 6 Gb
Memory to Store One Sparse Matrix: 1 Mb

NEVER AT ANY POINT should you use FULL MATRICES in the 
finite‐difference method.  Not even for intermediate steps.  NEVER!
Lecture 10 Slide 49

25