Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Numerical methods in Engineering
Typology: Lecture notes
1 / 20
This page cannot be seen from the preview
Don't miss anything!
Principle: derivatives in the partial differential equation are approximated
by linear combinations of function values at the grid points
u i ≈ (^) u ( x i ) ,
i = 0
grid points
x i = (^) i ∆ x
mesh size
NX
First-order derivatives
∂x ∂u
(^) (¯ x )
lim
∆ x → 0 u (¯ x (^) + ∆
x ) (^) −
(^) u (¯ x )
x
lim
∆ x → 0 u (¯ x ) (^) −
(^) u (¯ x (^) −
(^) ∆
x )
x
lim
∆ x → 0 u (¯ x (^) + ∆
x ) (^) −
(^) u (¯ x (^) −
(^) ∆
x )
x
(by definition)
Geometric interpretation
x i+
x i
x i 1
u
x
exa
t
en tral
forw
ard
ba
kw
ard
x
x
∂x∂u^
) i ≈
u i
− u i
∆ x
forward difference
∂x∂u^
) i ≈ u i − u
i− 1
∆ x
backward difference
∂x∂u^
) i ≈
u i
− u i − 1
2∆
x
central difference
Taylor series expansion
u ( x ) =
( x − x i ) n
n !
∂^ n u
∂x
n ) i (^) ,
u ∈
C ∞ ([
, X
1 :
u i
u i
(^) ∆
x ( ∂x∂u^
i
(∆
x ) 2
2 u
∂x
2 ) i
(∆
x ) 3
3 u
∂x
3 ) i
(^)...
2 :
u i − 1 =
u i (^) −
(^) ∆
x ( ∂x∂u^
i
(∆
x ) 2
2 u
∂x
2 )
i
−
(∆
x ) 3
3 u
∂x
3 )
i
(^)...
Central difference scheme
1
(^) T 2
2 u
∂x
2 ) i = u i
u i (^) +
(^) u i − 1
x ) 2
x ) 2
Alternative derivation
2 u
∂x
2 ) i = [ ∂
∂x
∂x∂u^
lim
∆ x → 0 ( ∂x∂u^
) i
/ 2 (^) −
(^) ( ∂x∂u^
) i − 1 / 2
x
u i+
− u i
∆ x
u i − u i − 1
∆ x
x
u i+
u i
(^) u i − 1
x ) 2
Variable coefficients
f (^) ( x ) =
d ( x ) ∂x^ ∂u
diffusive flux
∂x∂f^
i
f i
/ 2 − (^) f i − 1 / 2
x
d i+
/ 2 (^) u i+
− u i
∆ x
(^) d i − 1 / 2 (^) u i − u i − 1
∆ x
x
d i
/ 2 u i
d i
/ 2
(^) d i − 1 / 2 ) u i +
(^) d i − 1 / 2 u i
− 1
x ) 2
∂ 2 u
∂x∂y
∂x∂
∂y∂u^
∂y∂
∂x∂u^
)
„
∂ 2 u
∂x∂y
« i,j
=
“ ∂y∂u^
” i+
,j −
“ ∂y∂u^
” i− 1 ,j
2∆
x
(^) O
(∆
x ) 2
„ ∂y∂u^
« i
u i
,j
− (^) u i+
,j − 1
2∆
y
(^) O
(∆
y ) 2
„ ∂y∂u^
« i − 1 ,j =
u i − 1 ,j
(^) −
(^) u i − 1 ,j − 1
2∆
y
(^) O
(∆
y ) 2
x i+
x i (^1)
x i
x
y
y j 1
y j
y j +
Second-order difference approximation
2 u
∂x∂y
i,j
u i
,j
(^) u i
,j − 1 − (^) u i − 1 ,j
(^) u i − 1 ,j − 1
x ∆ y
x ) 2 , (^) (∆
y ) 2 ]
One-sided approximation
∂x∂u^
i ≈
αu
i (^) +
(^) βu
i
(^) γu
i
x
u i
(^) u i
x ( ∂x∂u^
i
(∆
x ) 2
2 u
∂x
2 ) i
(∆
x ) 3
3 u
∂x
3 ) i
(^)...
u i
(^) u i
x ( ∂x∂u^
i
(2∆
x ) 2
2 u
∂x
2 ) i
(2∆
x ) 3
3 u
∂x
3 ) i
(^)...
αu
i
βu
i
γu
i+
∆ x = α + β + γ
∆ x
u i
β (^) + 2
γ ) (^) ( ∂x∂u^
) i (^) +
∆ x 2 ( β
γ ) ( ∂^ 2 u
∂x
2 ) i
(^) O
(∆
x 2 )
Second-order accurate if
α (^) +
(^) β (^) +
(^) γ
= 0
β (^) + 2
γ = 1
β
γ = 0
α (^) =
2 3 (^) ,
β = 2
− 2 1
∂x∂u^
) i = − 3 u i
u i
− u i
2∆
x
x 2 )
One-sided approximation
2 u
∂x
2 ) i ≈
αu
i (^) +
(^) βu
i
(^) γu
i
x 2
u i
(^) u i
x ( ∂x∂u^
i
(∆
x ) 2
2 u
∂x
2 ) i
(∆
x ) 3
3 u
∂x
3 ) i
(^)...
u i
(^) u i
x ( ∂x∂u^
i
(2∆
x ) 2
2 u
∂x
2 ) i
(2∆
x ) 3
3 u
∂x
3 ) i
(^)...
αu
i
βu
i
γu
i
∆ x 2 = α + β + γ
∆ x 2
u i (^) +
β
γ
∆ x
( ∂x∂u^
) i
β
γ
2
∂^ 2 u
∂x 2 ) i
(^) O
(∆
x )
First-order accurate if
α (^) +
(^) β
(^) γ
= 0
β
γ = 0
β
γ = 2
α (^) = 1
, β = − 2 , γ
∂x∂u^
) i = u i − 2 u i
u i
∆ x 2
x )
Boundary value problem
2 u
∂x
f
in Ω = (
u (0) =
(^) u (1) = 0
One-dimensional mesh
x N
x 1
x 0
x i+
x i
x i (^1)
x N (^) 1
0
1
u i ≈ u ( x i ) , f i = f
x i ) x i = i ∆
x,
N 1
,
i = 0
Central difference approximation
x ) 2
(^) u i − 1 − 2 u i
u i+
(∆
x ) 2 = f i , ∀ i
u 0 = (^) u N
= 0
Dirichlet boundary conditions
Result: the original PDE is replaced by a linear system for nodal values
Linear system for the central difference scheme
i = 1
(^) u 0 − 2 u 1 + u 2 (∆
x ) 2
(^) f 1
i = 2
(^) u 1 − 2 u 2 + u 3 (∆
x ) 2
(^) f 2
i = 3
(^) u 2 − 2 u 3 + u 4 (∆
x ) 2
(^) f 3
(^) N
− (^1)
u N (^) − 2 − 2 u N (^) − 1
u N
(∆
x ) 2
(^) f N (^) − 1
Matrix form
Au
N (^) − 1 × N (^) − 1
u, F
N (^) − 1
x ) 2
u 1
u 2
u · 3
u N (^) − 1 (^) ,
f 1
f 2
f · 3
f N (^) − 1
The matrix
is tridiagonal and symmetric positive definite
invertible.
Non-homogeneous Dirichlet BC
u (0) =
g 0
only
changes
u 0 = g 0 ⇒ 2 u 1
(^) u 2
x ) 2 = f 1 + g 0
x ) 2
first equation
Non-homogeneous Neumann BC
∂x∂u
(^) (1) =
(^) g 1
only
changes
u N (^) +
(^) −
(^) u N (^) − 1
x = g 1 ⇒ u N
(^) +
u N (^) − 1
xg
1
(^) u N (^) − 1 − (^2) u N
(^) u N (^) +
x ) 2
(^) f N
u N (^) − 1
(^) u N
x ) 2
(^) f N
g 1
x
Non-homogeneous Robin BC
∂x∂u
(^) (1) +
(^) αu
(^) g 2
and
change
u N (^) +
(^) u N (^) − 1
x
(^) αu
N = g 2 ⇒ u N
(^) +
u N (^) − 1 (^) −
(^) 2∆
xαu
N
xg
2
(^) u N (^) − 1 (^) −
(^2) u N
(^) u N (^) +
x ) 2
(^) f N
u N (^) − 1 (^) + (1 +
(^) α ∆ x ) u N
x ) 2
(^) f N
g 2
x
Boundary value problem
(^) ∂ 2 u
∂x
2 −
∂ 2 u
∂y
f
in Ω = (
u = 0
on Γ =
Uniform mesh: ∆
x = ∆
h,
h 1
ix
ix 1
ix
h
h
x h h
1
1
0
y
ve-p oin t sten (^) il
j^ y 1 jy jy
u i,j
u ( x i , y
j (^) ) ,
f i,j
f (^) ( x i , y
j (^) ) ,
x i , y
j (^) ) = (
ih, jh
i, j
Central difference approximation
h 2 )
(^) u i − 1 ,j (^) + u i,j − 1 − 4 u i,j (^) + u i
,j (^) + u i,j
h 2
f i,j (^) ,
i, j
u i, 0
= (^) u
i,N
u 0 ,j
=
u N,j
i, j
2D Poisson equation
(^) ∂ 2 u
∂x 2 −
∂ 2 u
∂y 2 = (^) f
in Ω
g 0
on Γ
Difference equation − (^) u 1 (^) +
(^) u 2 − (^4) u 0 (^) +
(^) u 3
(^) u 4
h 2
f 0
curvilinear boundary
Ω
Q
P Γ
Ω
Q
0
4
h
2
1 3
stencil of^ R
(^) Q
δ Γ
Linear interpolation
u ( R ) =
u 4 ( h (^) −
(^) δ ) +
(^) u 0 δ
h = g 0 ( R ) ⇒ u 4 = − u 0 δ
h (^) −
(^) δ
(^) g 0 ( R )
h
h (^) −
(^) δ
Substitution yields
u 1 (^) +
(^) u 2 −
( 4 +
δ
h − δ ) u 0
(^) u 3 = h 2 f 0 +
(^) g 0 ( R )
h
h − δ
Neumann and Robin BC are even more difficult to implement
Purpose: to provide a simple treatment of curvilinear boundaries
dire
t (^) mapping
in v erse
(^) mapping
P
d
b
a
a
b
d
omputational
domain
artesian
grid
b o (^) dy-
tted
grid
ph
ysi
al
domain
nozzle
re
tangle
x
y
P
The original PDE must be rewritten in terms of (
ξ, η
) instead of (
x, y
) and
Derivative transformationsdiscretized in the computational domain rather than the physical one.
∂x ∂u
∂y ∂u
difficult to compute
∂ξ ∂u
∂η ∂u
easy to compute
Inverse mapping
x ( ξ, η
y ( ξ, η
Metrics transformations
∂x ∂ξ
∂y ∂ξ
∂x ∂η
∂y∂η
unknown
∂ξ ∂x
∂η ∂x
∂ξ ∂y
,
∂η∂y
known
Chain rule
∂ξ ∂u
∂x∂u
∂ξ∂x
∂y∂u
∂ξ∂y
∂η∂u
∂x∂u
∂η∂x
∂y∂u
∂η∂y
∂η∂u ∂ξ∂u
∂ξ∂x
∂ξ∂y
∂η∂x
∂η∂y
J
∂y∂u ∂x∂u
where
∂ ( x,y
)
∂ ( ξ,η
)
is the Jacobian which can be inverted using Cramer’s rule
Derivative transformations
∂x ∂u
det
∂ξ∂u^
∂η∂y
∂η∂u
∂ξ∂y
∂y∂u
det
∂η∂u^
∂ξ∂x
∂ξ∂u
∂η∂x
Total differentials for both coordinate systems
ξ ( x, y
η ( x, y
dξ
∂x∂ξ
(^) dx
∂y∂ξ (^) dy
dη
∂x∂η
(^) dx
∂y∂η (^) dy
dη dξ
∂x∂ξ
∂y∂ξ
∂x∂η
∂y∂η
dy dx
x ( ξ, η
y ( ξ, η
dx
∂ξ∂x (^) dξ
(^) +
∂η∂x (^) dη
dy
∂ξ∂y (^) dξ
∂η∂y (^) dη
dy dx
∂ξ∂x
∂η∂x
∂ξ∂y
∂η∂y
dη dξ
∂x∂ξ
∂y∂ξ
∂x∂η
∂y∂η
∂ξ∂x
∂η∂x
∂ξ∂y
∂η∂y
− 1 =
det
∂η∂y
∂η^ ∂x
∂ξ^ ∂y
∂ξ∂x
Relationship between the direct and inverse metrics
∂x ∂ξ
det
∂η ∂y
∂x∂η
det
∂ξ ∂y
∂y∂ξ
det
∂η ∂x
∂y∂η
det
∂ξ∂x
