Introduction to computational

fluid dynamics

Lecture 2 : CFD discretization

Navier Stokes Equation
⚫ Make the student understand the role of C in FD, its applicability,
potential and limitations
⚫ Give a basic foundation in numerical analysis, by teaching the
relevance of accuracy and stability
⚫ Give a working idea of the various choices of numerical methods
and discretization schemes by applying them to simple model
equations. In doing this, always remind them of the connection
with the big picture.
⚫ Make the student knowledgeable about the various terminologies
in practical CFD (Grids, BCs, Approximations, Schemes etc)
⚫ Ingrain the basics of good CFD practice (be aware of the
applicability/feasibility of a particular model, its limitations, choose
the right boundary conditions, ascertain grid/time independence,
verification/validation)
⚫ By the end of the class, the student should be in a position to set
up simple microfluidic CFD problems and analyze them
Navier-Stokes Equations
Conservation of Mass
 
+   V = 0
t

Conservation of Momentum

 u u u u  p    u 2      u v     u w 
  + u + v + w  = − +    2 −  V  +    +  +    +  + Bx
 t x y z  x x   x 3  y   y x  z   z x 
 v v v v  p    v 2      u v     v w 
  + u + v + w  = − +    2 −  V  +    +  +    +  + B y
 t x y z  y y   y 3  x   y x  z   z y 
 w w w w  p    w 2      w u     v w 
  + u +v + w  = − +    2 −   V  +    +  +    +  + Bz
 t x y z  z z   z 3   x   x z   y   z y 
Incompressible Navier-Stokes Equations
Conservation of Mass

 V = 0

Conservation of Momentum
 u   p
  + V   u  = − +  2 u + B x
 t  x
 v   p
  + V   v  = − +  2 v + B y
 t  y
 w   p
  + V   w  = − +  2 w + B z
 t  z

For isothermal flows, we have four unknowns: p, u, v, w.

Incompressible flow
Euler Equations
⚫ Neglecting all viscous terms in the Navier-Stokes
equations yields the Euler equations:
 
+   V = 0
t
u  p
+   Vu = − + Bx
t x
v  p
+   Vv = − + B y
t y
w  p
+   Vw = − + Bz
t z
⚫ No transport properties (viscosity or thermal conductivity) are needed.
⚫ Momentum and energy equations are greatly simplified.
⚫ But we still have five unknowns: , p, u, v, w.
⚫ The Euler equations provide a reasonable model of compressible fluid flows at
high speeds (where viscous effects are confined to narrow zones near wall
boundaries).
Incompressible NS
Exact solutions

 u   p
  + V  u  = − +  2u + Bx
 t  x

⚫ Incompressible viscous steady flow between two

parallel plate
d 2u dp
 2 =
dy dx
Integrating twice, we obtain

⚫ u=
1 dp 2
2 dy
(
y − a2 )
From Mathematical Models to Numerical
Simulations
Sphere Motion in Fluid Flow
Sphere Motion in Fluid Flow (Matlab)
clear
rho=1000;
Cd=1;
m=5;
r=0.05;
fac=rho*Cd*pi*r^2/(2*m);
v=1;
save dragpar ;
x=[0:0.1:10];
%step size
h=1.0;
% Euler's method, forward finite difference
t=[0:h:10];
N=length(t);
u_e=zeros(N,1);
x_e=zeros(N,1);
u_e(1)=0;
x_e(1)=0;
for n=2:N
u_e(n)=u_e(n-1)+h*fac*(v^2-2*v*u_e(n-1)+u_e(n-1)^2);
x_e(n)=x_e(n-1)+h*u_e(n-1);
end
% Runge Kutta
u0=[0 0]';
[tt,u]=ode45(@dudt,t,u0);
Partial Differential Equations
Linear in two variables:
 2u  2u  2u
A( x, y ) 2 + B( x, y ) + C ( x, y ) 2 + ...... = 0
x xy y
Usual classification at a given point (x,y):
B 2 − 4 AC  0 → elliptic Laplace' s equation
B 2 − 4 AC = 0 → parabolic Heat conduction equation
B 2 − 4 AC  0 → hyperbolic Wave equation
From the numerical point of view
Initial Value Problem ( time evolution) Hyperbolic or Parabolic
Boundary Value Problem ( static solution) Elliptic

Computational Concern:
Initial Value Problem : Stability
Boundary Value Problem: Efficiency

13
Nondimensionalization of the NSE
⚫ Purpose: Order-of-magnitude analysis of the terms in
the NSE, which is necessary for simplification and
approximate solutions.
⚫ We begin with the incompressible NSE

⚫ Each term is dimensional, and each variable or

property (, V, t, , etc.) is also dimensional.
⚫ What are the primary dimensions of each term in the
NSE equation?
Nondimensionalization of the NSE
⚫ To nondimensionalize, we choose scaling
parameters as follows
Nondimensionalization of the NSE
⚫ Next, we define nondimensional variables, using the
scaling parameters in Table 10-1

⚫ To plug the nondimensional variables into the NSE, we

need to first rearrange the equations in terms of the
dimensional variables
Nondimensionalization of the NSE
⚫ Now we substitute into the NSE to obtain

⚫ Every additive term has primary dimensions

{m1L-2t-2}. To nondimensionalize, we multiply
every term by L/(V2), which has primary
dimensions {m-1L2t2}, so that the dimensions
cancel. After rearrangement,
Nondimensionalization of the
NSE
⚫ Terms in [ ] are nondimensional parameters

Strouhal number Euler number Inverse of Froude Inverse of Reynolds

number squared number

Navier-Stokes equation in nondimensional form

Nondimensionalization of the NSE
⚫ Nondimensionalization vs. Normalization
⚫ NSE are now nondimensional, but not necessarily normalized.
What is the difference?
⚫ Nondimensionalization concerns only the dimensions of the
equation - we can use any value of scaling parameters L, V, etc.
⚫ Normalization is more restrictive than nondimensionalization.
To normalize the equation, we must choose scaling parameters
L,V, etc. that are appropriate for the flow being analyzed, such
that all nondimensional variables are of order of magnitude
unity, i.e., their minimum and maximum values are close to 1.0.

If we have properly normalized the NSE, we can compare the relative

importance of the terms in the equation by comparing the relative
magnitudes of the nondimensional parameters St, Eu, Fr, and Re.
Creeping Flow
⚫ Also known as “Stokes Flow” or “Low
Reynolds number flow”
⚫ Occurs when Re << 1
⚫ , V, or L are very small, e.g., micro-
organisms, MEMS, nano-tech, particles,
bubbles
⚫  is very large, e.g., honey, lava
Inviscid Regions of Flow
⚫ Definition: Regions where net viscous forces
are negligible compared to pressure and/or
inertia forces
~0 if Re large

Euler
Equation
Finite Differences Method
f
f i +1 = f ( xi +1 ) = f ( xi + h) = f ( xi ) + f ( xi ) h + f ( xi ) h / 2!+....... 1
2
fi-1 fi fi+1
solve for first derivative
f ( xi + h) − f ( xi )
f ( xi ) = − f ( xi ) h / 2!+......(error of order h )
h x
xi-1=x-h xi xi+1=x+h
f −f
f ( xi ) = i +1 i (forward difference )
h

f i −1 = f ( xi −1 ) = f ( xi − h) = f ( xi ) − f ( xi ) h + f ( xi ) h 2 / 2!+....... 2
solve for first derivative
f ( xi ) − f ( xi − h) combining formulas 1 and 2
f ( xi ) = + f ( xi ) h / 2!+......(error of order h )
h we obtain the central
difference approximation
f i − f i −1
f ( xi ) = (backward difference ) which is second order
h
f i +1 − f i −1
f ( xi ) =
2h

Boundary Value P 22
Second derivative formula
f i +1 = f i + f i h + f i h 2 / 2!+ f i h3 / 3!+......... 1

f i −1 = f i − f i h + f i h 2 / 2! − f i h3 / 3! + ......... 2

Add the two equations

f i +1 + f i −1 = 2 f i + f i h 2 + (h 4 ) and solving for the second derivative

f i +1 − 2 f i + f i −1
f ( xi ) = (central difference )
h2

Boundary Value P 23
Derivative formulas
f i +1 − f i −1
f ( xi ) =
2h
f i +1 − 2 f i + f i −1

f ( xi ) =
h2 Central Difference Expressions
f − 2 f i +1 + 2 f i −1 − f i − 2 with error of order h2
f ( xi ) = i + 2
2h 3
f − 4 f i +1 + 6 f i − 4 f i −1 + f i − 2
f ( xi ) = i + 2
h4
− f i + 2 + 4 f i +1 − 3 f i
f ( xi ) =
2h
− f i +3 + 4 f i + 2 − 5 f i +1 + 2 f i
f ( xi ) =
h2
Forward Difference
− 3 f i + 4 + 14 f i +3 − 24 f i + 2 + 18 f i +1 − 5 f i
f ( xi ) = Expressions
2h 3 with error of order h2
− 2 f i +5 + 11 f i + 4 − 24 f i +3 + 26 f i + 2 + 3 f i +1 + 3 f i
f ( xi ) =
h4
Boundary Value P 24
Derivative formulas 2

3 f i − 4 f i −1 + f i − 2
f ( xi ) =
2h
2 f − 5 f i −1 + 4 f i − 2 − f i −3
f ( xi ) = i
h2
5 f − 18 f i −1 + 24 f i − 2 − 14 f i −3 + 3 f i − 4 Backward
f ( xi ) = i
2h 3 Difference
3 f − 14 f i −1 + 26 f i − 2 − 24 f i −3 + 11 f i − 4 − 2 f i −5 Expressions
f ( xi ) = i
h4 with error of order
h2

Boundary Value P 25
Finite Differences Stencil
S

R u xx + u yy = 0 or ( u = 0 ) on R

Using Central Finite Difference approximation of

order h2 over a uniform grid of any size h we obtain the
following stencil to be applied at all internal points
1

1 -4 1

PDEs & Elliptic problems 26

Biharmonic equation Stencil
S
D   u = pz on R
R Here D is the stiffness of a plate, pz is the load in
a direction perpendicular to the plate, and the
= 2+ 2
 
2 2
operator
x y
1

2 -8 2

 pz 
1 -8 20 -8 1 =h  
4

 D  m,n
2 -8 2
Pivotal point m,n

1
PDEs & Elliptic problems 27
Laplace Equation: Unequal grid
u2 ( x , y+h2 )
( x-h3, y ) h3 h2 ( x+h1, y ) We
approximate Laplace operator by:
u3 u0 ( x , y ) h1 u1
uxx + uyy = a0u0 + a1u1 + a2u2 + a3u3 + a4u4 ( I )
h4 In order to determine the coefficients ai we use
u4 the expansions:
( x , y-h4 )

h12 h32
u1 = u0 + h1u x + u xx + (h13 ) u3 = u0 − h3u x + u xx + (h33 )
2 2
h22 h42
u2 = u0 + h2u y + u yy + (h23 ) u 4 = u0 − h4u y + u yy + (h43 )
2 2
Substituting these expansions in ( I ) produces:
u xx + u yy = u0 (a0 + a1 + a2 + a3 + a4 ) + u x (a1h1 − a3h3 ) + u y (a2 h2 − a4 h4 ) +
u xx (a1h1 + a3h3 ) + 12 u yy (a2 h2 + a4 h4 ) + third order error terms
1 2 2 2 2
2

Validity of this equation requires the satisfaction of a set of equations shown next →

PDEs & Elliptic problems 28

Laplace Equation: Unequal grid 2
Validity of the 1 1 1 1 1  a0   0 
representation requires the     
satisfaction of the set of  0 h1 0 − h3 0  a1   0 
equations shown → 0 0 h2 0 − h4  a2  =  0 
    
 0 h12 0 h32 0  a3   2 
0 0 h22 h42  a4   2 
 0
The solution of this system is:
 1 1  2 2
a0 = −2  +  , a = , a = ,
h1 (h1 + h3 ) h2 (h2 + h4 )
1 2
 h1h3 h2 h4 
2 2
a3 = , a4 =
h3 (h1 + h3 ) h4 (h2 + h4 )

Note that for equal h’s this reduces to:

1
2
(− 4u0 + u1 + u2 + u3 + u4 )
h
PDEs & Elliptic problems 29
Example
15
(2,5) u is determined on the boundaries by the function:
4 u = x2 – y2 .
(0,4) 14
3 Obtain u at all points so that it satisfies Laplace’s equation.
13
1 2 12 In determining the coordinates of the points on the side line,
5 please note that the equation of the line is: y = -(5/5 )(x – 7)
11

6 10
7 8 9 (7,0)

A) First number all points, starting with the interior points

B) Then evaluate u at the boundary points

PDEs & Elliptic problems 30

Example Cont. 1 Values of u on the boundary:
(2,5)
15 u4 = -16
4 u5 = -4 u11 = 36 – 1 = 35
14
(0,4) 3
13 u6 = 0 u12 = 25 – 4 = 21
5 1 2 12 u7 = 4 u13 = 16 – 9 = 7
11 u8 = 16 u14 = 9 – 16 = -7
6 10 u9 = 36 u15 = 4 – 25 = -21
7 8 9 (7,0)
u10 = 49
c) Write the equations at the internal points:
− 4u1 + u2 + u3 + (−4) + 4 = 0 → − 4u1 + u2 + u3 = 0
1 1 1 1 1 1 1
−  + u2 +  21 +  7 + u1 + 16 = 0 → − u2 + u1 = −12
2 2 3 3 6 6 6
1 1 1 1 1 1 1
−  + u3 +  − 7 +  −21 +  −16 + u1 = 0 → − u3 + u1 = 12
2 2 3 3 6 6 6

PDEs & Elliptic problems 31

Example Cont. 2 Values of u on the boundary:
(2,5)
15 u4 = -16
4 u5 = -4 u11 = 36 – 1 = 35
14
(0,4) 3
13 u6 = 0 u12 = 25 – 4 = 21
5 1 2 12 u7 = 4 u13 = 16 – 9 = 7
11 u8 = 16 u14 = 9 – 16 = -7
6 10 u9 = 36 u15 = 4 – 25 = -21
7 8 9 (7,0)
u10 = 49
 − 4 1 1  u1   0 
The equations of the previous  1    
slide can be written as:  6 − 1 0  u2  =  − 12 
 1 −  u   12 
 6 0 1 3   
Direct solution of these equations produces:
u1 = 0, u2 = 12, u3 = -12
Note that the system is diagonally dominant

PDEs & Elliptic problems 32

Parabolic PDEs
area A
insulated rod
heat flux in heat flux out

T(x,t) T(x+x,t)
time line
x

T = 0o
What is the T = 0o T
heat flux in = − k A
at all times temperature at all times
x x ,t

later ? T
heat flux out = − k A
x=0 T = T0 f ( x ) x=L x axis x x + x ,t
at time = 0
Balance equation:
T T
−
T T T  T k x x
cv A x = −k A +k A = x + x ,t x ,t
which can be written as :
t x x ,t x x + x ,t  t cv x

in the limit
 T k  2T
→ =
 t cv x 2
PDEs & Parabolic problems 33
nondimensionalization
 T k  2T
= select nondimensional variables x, t, T using
 t cv x 2
x t T
x = , t= , T=
L t scale T0
and the equation becomes:
T0  T k T0  2T
=
t scale  t cv L2 x 2
cv L2  T  2T
choosing t scale = we obtain: = 2
k  t x

PDEs & Parabolic problems 34

 Numerical Solution
consider a grid:
u = g1(t) u = g2(t)
u t = u xx
 
t = k
 x=h
u ( x − h, t ) − 2u ( x, t ) + u ( x + h, t )
x=0 u = f(x) x=L
u xx ( x, t ) =
h2
u ( x, t + k ) − u ( x, t )
ut ( x, t ) =
k
this results in an explicit method:
k
2
u (x − h, t ) − 2u (x, t ) + u (x + h, t )
u ( x, t + k ) = u ( x, t ) +
h
or u ( x, t + k ) =  u (x − h, t ) + (1 − 2 ) u ( x, t ) +  u ( x + h, t )
k
=
where h2

PDEs & Parabolic problems 35

Stencil for explicit solution

2
k
k = 2
h
3 0 h 1
h

u2 =  u1 + (1 − 2 ) u0 +  u3

PDEs & Parabolic problems 36

Stable solution- the Implicit method

(x-h , t) (x , t) (x+h , t)

3 h 0 h 1
The differenti al equation ut = u xx translates to
k
(x , t - k) u0 − u4 u3 − 2u0 + u1
4 = 2
or
k h
u0 − u4 =  (u3 − 2u0 + u1 ) where  = k 2
h

this produces the diagonally dominant

tridiagona l system at each row of points :
λ u3 − (1 + 2 λ )u0 +  u1 = −u4

PDEs & Parabolic problems 37

Example on Implicit method
t
7 8 k=h=1/3 , =3 , 1 + 2 = 7
5 6

u=0 3 4 u=1 for the points 1 and 2 at the first level

up we obtain
1 2

x
3 * 0 - 7 u1 + 3 u2 = - 0.3
u =.3 u =.7
3 u1 - 7 u2 + 3 * 1 = - 0.7

using one decimal digit arithmetic the

solution of this system is:
u2 = .7 and u1 = .3

Stable and consistent method.

PDEs & Parabolic problems 38

Crank - Nicholson Method
u u0 − u 4
3 0 1 = central difference
t A k
A k
 2u u3 − 2u0 + u1
=
7 4 8 x 2 0 h2
h
 2u u7 − 2u4 + u8
=
x 2 4 h2
 2u 1   2u  2u  u3 − 2u0 + u1 + u7 − 2u4 + u8
=  2 + 2 =
x 2 A
2  x 0 x 
4
2 h 2

u0 − u 4 =
k
2
u3 − 2u0 + u1 + u7 − 2u4 + u8  using =
k
2h h2

 u3 − 2( + 1) u0 +  u1 = −2u4 −  (u7 − 2u4 + u8 )

This is more accurate without
increasing the work considerably.
Same form LHS of the equation

PDEs & Parabolic problems 39

 A Boundary Value Technique 1
t t
u = U (x)
t =T

u = g1 (t ) u = g 2 (t ) u = g1 (t ) u = g 2 (t )

u xx = ut + F ( x, t , u ) u xx = ut + F ( x, t , u )

x x
x=0 x=a x=0 x=a
u = f (x) u = f (x)

here U(x) is the

solution of the U xx = lim F ( x, t , U )
t →
problem for t → U (0) = lim g1 (t )
t →
large
U (a ) = lim g 2 (t ) solve the problem
summarized on the t →
and get U(x)
right
PDEs & Parabolic problems 40
A Boundary Value Technique 2

we can then write the finite difference analog of the equation

using the stencil

2 u(x,t+k)

u(x-h,t) 3 0 u(x,t) 1 u(x+h,t)

4 u(x,t-k)

see example next slide

PDEs & Parabolic problems 41

B V T Example 1
u xx − ut = x u x 0  x 1
the exact solution of this problem is:
u (0, t ) = 0 t0 u ( x, t ) = x e − t
u (1, t ) = e −t t0
( differentiate and substitute
u ( x,0) = x 0  x 1
to prove that this is indeed
the case)
As time tends to large values the problem simplifies to:

Uxx = x Ux with U(0) = 0 and U(1) = 0

which has the trivial solution U = 0 everywhere

We choose T = 2 ( this is our approximation to large time)

and we use k = h = 1/3

PDEs & Parabolic problems 42

B V T Example 2
t u=0 u=0

u=0
(1/3,2) (2/3,2)
(0,2) (1,2) u=e -2 =0.14

u=0
(0,5/3) (1,5/3) u=e -5/3
we now write the finite difference
9 10 equations for each of the 10
u=0
points and we obtain 10 equations
(0,4/3) to be solved simultaneously. This
(1,4/3) u=e -4/3
7 8 produces a solution which agrees
with the exact one in two decimal
u=0 places
(0,1) (1,1) u=e -1
5 6
The first two equations
u=0 ( points 1 and 2 ) are:
(0,2/3) (1,2/3) u=e -2/3
3 4
17 3 1
u=0
− 18u1 + u 2 − u3 = −
(0,1/3) (1,1/3) u=e -1/3
2 2 2
1 2
3 − 13
− 18u2 − u4 + 10u1 = −1 − 8e
u=0
u=1 x 2
(0,0)
(1/3,0) (2/3,0) (1,0) verify this
u=1/3 u=2/3

PDEs & Parabolic problems 43

Two Dimensional Diffusion Equation
u  2u  2u
= 2+ 2 Try to implement Crank-
t x y Nicholson Scheme on this
6problem
n+1 2 n+1 5 n+1

3n+1 0n+1
1n+1 x=y=h
7n+1 4n+1
A
8n+1 6n 2n 5n t = k

3n 0n 1n

7n 4n 8n

1 2 u1n+1 − 2u0n+1 + u3n+1 1 2 u 2n+1 − 2u0n+1 + u4 n+1

 u =
2 x 0 n+1 2
, 2
 y u0n+1 =
h h h h2
u0n+1 − u0n
k
1
(
= 2  x2 u0n+1 +  x2 u0n +  y2 u0n+1 +  y2 u0 n
2h
)
We must write one such equation per each interior point and then solve the
system at each time step.
PDEs & Parabolic problems 44
Two Dimensional Diff. Eqn Cont.

u  2u  2u
= 2+ 2
t x y u0n+1 − u0n
k
=
1
2h 2
(2
x u 0 n+1 +  x u n +  y u n +1 +  y u n
2
0
2
0
2
0
)
We must write one such equation per each interior point and then solve the
system at each time step.

Difficulty: Large system of equations which is sparse but not tridiagonal or even banded

Remedy: ADI method ( Alternating Direction Implicit solution)

a) write the above equation in the form:
u nj ,+l 1 = u nj ,l +
k
2h 2
(2 n +1
x u j ,l +  2 n
x u j ,l +  2 n +1
y u j ,l +  2 n
y u j ,l )
b) divide each time step into two steps of size t/2. In each sub step,
treat a different dimension implicitly i.e.

u nj ,+l 2 = u nj ,l +
1 k
2h 2
 (
2 n + 12
x u j , l +  2 n
y u j ,l ) and u nj ,+l 1 = u nj ,+l 2 +
1 k
2h 2
(
 2 n + 12
x u j , l +  2 n +1
y u j ,l )
Advantage: at each sub step only a tridiagonal system is solved
PDEs & Parabolic problems 45
Boundary Value Problems – an explicit method & stability

3k 9 10 11 12
ut t = u x x
u ( x,0) = f1 ( x) 
2k 5 6 7 8  I .C.
ut ( x,0) = f 2 ( x)
k 1 2 3 4
u (0, t ) = g1 (t ) 
 B.C
u (1, t ) = g 2 ( x)
h 2h 3h 4h 5h x
for the second row use the stencil:
for the first row use
2 ( x, t+k)
u ( x, k ) − u ( x,0)
= f 2 ( x) to obtain : ( x-h, t) 3 0 ( x, t) 1 ( x+h, t)
k
ui = f1 (i h) + k f 2 (i h) for i = 1,2,3,4
4 ( x, t-k)

u2 − 2u0 + u4 u3 − 2u0 + u1 k2
to obtain 2
= 2
or u2 = 2u0 − u4 + 2 u3 − 2u0 + u1 
k h h
 BVP – Example of explicit m.
11 12 13 14 15 ut t = u x x h= 1
6 , k= 1
2

6 7 8 9 10 u ( x,0) = x 
u=0 u=1
 I .C.
1 2 3 4 5
ut ( x,0) = 1
u (0, t ) = 0
 B.C
h 2h 3h 4h 5h 1 x
u (1, t ) = 1 
u=x , ut=1
From the formula ui = f1 (i h) + k f 2 (i h) we obtain:

u1 = 1
6 + 12 = 2 3
u2 = 1
3 + 12 = 5 6
u3 = 1
2 + 12 = 1
u4 = 2 3 + 1 2 = 7 6
u5 = 5 6 + 1 2 = 4 3
PDE's and the Wave
Equation 47
BVP – Ex of explicit method 2
ut t = u x x h= 1
6 , k= 1
2
11 12 13 14 15
u ( x,0) = x 
 I .C.
6 7 8 9 10 ut ( x,0) = 1
u=0 u=1
u (0, t ) = 0
 B.C
1 2 3 4 5 u (1, t ) = 1  2 ( x, t+k)

( x-h, t) 3 0 ( x, t) 1 ( x+h, t)
h 2h 3h 4h 5h 1 x
u=x , ut=1 4 ( x, t-k)

k2
From the formula u 2 = 2u0 − u 4 + 2 u3 − 2u0 + u1  ( stencil numbering) we obtain the
h
second row values at points 7,8,9 without using the boundary conditions

5 1 2 5  4
u7 = 2 • − + 9 − 2 • + 1 = then in the third row we can obtain:
6 3 3 6  3
3 4 3 5
3 u13 = 2 • − 1 + 9  − 2 • +  = 2
u8 = 2 3 2 3
2
It is very unreasonable to obtain u13 without using the boundary
5
u9 = conditions ( Why???) → Stability condition: k <= h
3
PDE's and the Wave
Equation
Boundary Value Problems – a stable implicit method
a) Find the values of u at the first row using the way described in the
explicit method earlier.

b) For each internal point in the first row write the FDE of the wave equation as
follows:
 2u u6 − 2u2 + u5 
= 
x 22
h2   2u 1   2u  2u 
  =  2 + 2 
6 2 5
u2
u7 − 2u4 + u8  x 0 2  x 2 x
2
4
= 
3 0 1 x 2 4 h2 
7 4 8
and
 2u u2 − 2u0 + u4
=
t 2 0 k2
this produces a tri diagonal system of equations for the values of u at the second row

 h2   h2  h2
u6 − 21 + 2 u2 + u5 = - u7 + 21 + 2 u4 − u8 − 4 2 u0
PDE's and the Wave
 k   k  k
Equation
Hyperbolic PDE: B2 -4 A C > 0
Partial Differential Equations Hyperbolic PDE
Partial Differential Equations Hyperbolic PDE
Partial Differential Equations Hyperbolic PDE
Partial Differential Equations Hyperbolic PDE
Hyperbolic PDE Method of Characteristics
Hyperbolic PDE
Start of Integration: Euler and Higher Order starts
Matlab coding
Courant-Fredrichs-Lewy Condition (1920’s)

⚫ Basic idea: the solution of the Finite-Difference (FD)

equation can not be independent of the (past) information
that determines the solution of the corresponding PDE
⚫ In other words: “Numerical domain of dependence of FD
scheme must include the mathematical domain of
dependence of the corresponding PDE”

CFL NOT satisfied CFL satisfied

CFL: Linear convection (Sommerfeld Eqn)
Example
CFL: 2nd order Wave equation Example
Elliptic PDE: Poisson Equation
Elliptic PDE: Poisson Equation
Matlab code