Download Numerical Solution of Differential Equations using Finite Difference Method and more Slides Engineering Analysis in PDF only on Docsity!
Numerical Solutions of theNumerical Solutions of the
Diffusion EquationDiffusion Equation
Larry Caretto
Mechanical Engineering 501AB
Seminar in Engineering Analysis
March 18, 2009
2
Outline
- Review last class
- Numerical solutions of the diffusion
equation in one space dimension
- Explicit algorithm
- Stability of algorithms
- Crank-Nicholson algorithm
- (Fully) implicit algorithm
- DuFort-Frankel algorithm
3
Review Numerical Analysis
- Transform differential equation into a
system of algebraic equations
- Obtain solution for discrete points in
domain
- Two basic approaches: finite differences
and finite elements
- Start with finite elements
- Get expressions for derivatives and
measure of error with their use
4
Review Finite Difference Grids
- Grid notation for four independent
variables: x, y, z and t
x
= x
min
x
N
= x
max
x
i
i-
= Δx
i
y
= y
mjn
y
M
= y
max
y
j
j-
= Δy
j
z
= z
mkn
z
K
= z
max
z
k
k-
= Δz
k
t
= t
min
t
L
= t
max
t
n
n-
= Δt
n
- Dependent variable u(x,y,z,t) at discrete
points u(x
i
, y
j
, z
k
, t
n
)
- Use notation below for this value of u
( , , , ) i j k n
n
ijk
u = u x y z t
5
Review Truncation Error
- If we truncate series after m terms
∑ ∑
∞
= (^) = = + =
= + 0 1
( -) !
1 (-) !
1 ( ) nm
n
xa
n
m n
n
n
xa
n
n
xa dx
d f
n
xa dx
d f
n
fx
Terms used Truncation error, ε
m
- Can write truncation error as single term
at unknown location (derivation based
on the theorem of the mean)
1 1
1
1
=
∞ +
= (^) ∑ =
m
x
m
m
nm
n
xa
n
n
m x a
dx
d f
m
x a
dx
d f
n
ξ
ε
6
Review Order of the Error
- Derivative expressions have error that is
proportional to h
n
- This power, n, is called the order of the
error
n
) to indicate this error
- Reducing step size by a factor of a
reduces n
th
order error by a
n
n
h
h
⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛ ≈
1
ε ε
7
Review Derivative Expressions
- First-order error, first derivatives
( ) 2
' 1 1 2 Oh h
f f f
i i i +
−
' (^1)
Oh
h
f f
f
i i i
h
f f f
i i i +
−
−
- Second-order error, first derivatives
Oh
h
f f f
f
i i i
i
Oh
h
f f f
f
i i i
i
'' 1 1 2 Oh
h
f f f f
i i i i +
8
Find f’ and f’’ for sin at x = 1
( ) 2
' 1 1 2 Oh h
f f f
i i i +
−
Second order central
( )
2 2
'' 1 1 2 Oh h
f f f f
i i i i +
2 (. 1 )
' sin(^1 +.^1 )−sin(^1 −.^1 ) f (^) i =^2 (. 1 )
sin( 1. 1 ) sin( 1. 1 ) 2 sin( 1 ) ''
f
2 (. 01 )
' sin(^1 +.^01 )−sin(^1 −.^01 )
i
f
2 (. 001 )
' sin(^1 +.^001 )−sin(^1 −.^001 ) fi =
2 (. 01 )
sin( 1. 01 ) sin( 1. 01 ) 2 sin( 1 ) ''
2 (. 001 )
sin( 1. 001 ) sin( 1. 001 ) 2 sin( 1 ) ' '
9
Review Roundoff Error
- Possible in derivative expressions from
subtracting close differences
x
: f’(x) ≈ (e
x+h
x-h
)/(2h)
and error at x = 1 is (e
1+h
1-h
)/(2h) - e
3
718282 4. 510 2 ( 0. 1 )
004166 2. (^722815) − − =
− E = x
9
718281828459 4. 510 2 ( 0. 0001 )
7185536702 2. (^7180100139) − − =
− E = x
9
718281828 5. 910 2 ( 0. 0000001 )
71828210028724 2. (^71828155660388) − − =
− E = x
Second order error
Figure 2-1. Effect of Step Size on Error
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E+
1.E+
1.E-17 1.E-15 1.E-13 1.E-11 1.E-09 1.E-07 1.E-05 1.E-03 1.E-
Step Size
Error
Slope of a log error – log h
plot is order n
Δ(log ε) = 10; Δ(log h) = 5
Slope = order = 2
11
Numerical PDE Solutions
- Define a finite-difference grid in the
independent variables (x, y, z, t)
- Place grid points on region boundary
whose values are found from boundary
conditions for the problem
- At some grid location convert differential
equation into a finite difference equation
- Observe truncation error in process
- Neglect truncation error to get set of
algebraic equations to solve
12
Diffusion Equation
- Apply difference formulas derived for
ordinary derivatives to partial derivatives
- Use notation to consider different
coordinate directions
- Apply to diffusion equation
- Grids xi = x 0 + iΔx and tn = t 0 + nΔt
- Try finite difference expressions below
to get simple finite-difference equation
n
i
x
u
t
u
α
[( ) ] ( )
2 ( )
2 2
1 1 2
1 2
O x x
u u u
x
u O t and t
u u
t
u
n i
n i
n i
n
i
n i
n i
n
i
∂
−
∂
∂ (^) + −
19
Explicit Results f = 0.
Error t = 0.20 0 385.9 422.5 385.9 0
Exact t = 0.20 0 125.1 176.9 125.1 0
n = 5 t = 0.20 0 -260.9 599.3 -260.9 0
n = 4 t = 0.16 0 419.2 -224.2 419.2 0
n = 3 t = 0.12 0 -35.3 639.6 -35.3 0
n = 2 t = 0.08 0 539.2 180.8 539.2 0
n = 1 t = 0.04 0 360 1000 360 0
n = 0 t = 0+ 0 1000 1000 1000 0
t = 0 1000 1000 1000 1000 1000
x = 0.00 x = 0.25x = 0.50x = 0.75 x = 1.
i = 0 i = 1 i = 2 i = 3 i = 4
20
What Happened?
- We are seeing effects of instability
- Difference equations may not converge
- Unstable equations grow without bound
- May have stable equations that produce
incorrect results
- Conditional stability requires step size less
than that needed for accuracy
- Goal of absolute stability not always possible
- Discussions of stability complex, can
sometimes use physical arguments
21
Stability of Explicit Method
i+
and u
i-
are fixed an
increase in u
i
n
should increase u
i
n+
- If f is greater than 0.5, an increase in u
i
n
will cause a decrease in u (^) i
n+
- We can avoid this incorrect result by
keeping f = αΔt/(Δx)
≤ 0.
- This imposes a time step limit that may
be less than the limit required for
accuracy in the solution
22
FTCS Truncation Error
- Derivation in appendix for notes on
solving PDEs gives this equation
∞
=
− −
2
2
1 2 2 2
k
n
i
k
k k n k i
x
T
k
f
k
TE α x
+ L
n
i
n
i
n
i
x
T
f
x
x
T
f
x
TE
= 1/6 eliminates
first term in the truncation error
23
Crank-Nicholson Method
- Seek more accurate time derivative
- Provides implicit method
n+
depends on other u
n+
- More work per step, but can take longer
time steps with this method
- Apply to diffusion equation at time n + 1/
[( )] [( )]
n
i
n
i
n
i
n
i
n
i
n
i x
u
O t
t
u u
O t
t
u u
t
u
24
Space Derivative at t
n+1/
- Take average of space derivative at
time steps n and n + 1
- Show average is second order accurate
( ) 2
..... 2 4 48
..... 24
2 2
2 2
..... 2 6
..... 2 6
1 1 2
4 ''''
2 1 1 ''
4 ''''
2 ''' 1 1
3 '''
2 ' '' 1
3 '''
2 ' '' 1
O h
h f f f
h f
f f f
h f
h f f f f
h f
h f f fh f
h f
h f f fh f
i i i i
i i i
i i i i i
i i i i i
i i i i i
− − + =
=
= − + − +
= + + + +
−
+
25
Using Space Derivative at t
n+1/
• Apply average to space derivative
[( ) ]
O t
x
u
x
u
x
u
n
i
n
i
n
i
[( ),( )] 0 ( )
2
( )
2
2
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
Δ
Δ
α + − −
Δ
−
∂
∂ −α ∂
∂
O t x x
u u u
x
u u u
t
u u
x
u
t
u
n i
n i
n i
n i
n i
n i
n i
n i
n
i
n
i
• Substitute into diffusion equation
• Introduce f = αΔt/(Δx)
and rearrange
26
Crank-Nicholson Equation
• Resulting equation has three values at
new time step
[ ]
n
i
n
i
n
i
n
i
n
i
n
i u u fu
f
u
f
u fu
f
• Tridiagonal system of equations easily
solved by special application of Gauss
elimination called Thomas algorithm
n
i
n
i
n
i
n
i
− fu + + fu − fu = R
[ ]
n
i
n
i
n
i
n
i
n
i
n
fui 2 ( 1 f ) u fu fu 1 u 1 2 ( 1 f ) u
27
Crank-Nicholson Equations
−
−
−
−
n N
n N
n N
n
n
n n
n N
n N
n
n
n
R fu
R
R
R
R fu
u
u
u
u
u
f f
f f
f f
f f f
f f f
f f
1
2
3
2
1 0
1 1
1 2
1 3
1 2
1 1
M
M
M
M
L
L
M M M M O M M
L
L
L
• Consider case where boundary
potentials u
and u
N
are specified
• Rewrite equations in matrix form to
show tridiagonal structure
28
Thomas Algorithm
• General format for tridiagonal equations
N
N
N
N
N N
N N
D
D
D
D
D
x
x
x
x
x
A B
B C
A B
A B C
A B C
B C
M
M
M
M
L
L
M M M M O M M
L
L
L
29
Thomas Algorithm II
N
N
N
N N
F
F
F
F
F
x
x
x
x
x
E
E
E
E
M
M
M
M
L
L
M M M M O M M
L
L
L
• Gauss elimination upper triangular form
30
Thomas Algorithm III
• Forward computations
- Initial: E 0 = – C 0 / B 0 F 0 = D 0 / B (^0)
- For i = 1,… N-1:
i i i
i i i i
i i i
i i
B AE
D AF
F
B AE
C
E
N N N
N N N
N N
B AE
D AF
x F
• Back substitute: x
i
= F
i
+ E
i
x
i+
• Get last x
value first
37
Fully Implicit Results III
Error t = 0.0125 0 0.779 1.542 2.273 2.
Exact t = 0.0125 0 50.43 100.66 150.48 199.
n = 25 t = 0.0125 0 51.21 102.20 152.76 202.
n = 24 t = 0.012 0 52.30 104.36 155.96 206.
n = 23 t = 0.0115 0 53.46 106.67 159.38 211.
n = 22 t = 0.011 0 54.70 109.13 163.04 216.
n = 21 t = 0.0105 0 56.03 111.78 166.95 221.
n = 20 t = 0.01 0 57.47 114.62 171.16 226.
n = 19 t = 0.0095 0 59.02 117.70 175.71 232.
n = 18 t = 0.009 0 60.70 121.03 180.64 239.
x = 0 x = .01 x = .02 x = .03 x =.
i = 0 i = 1 i = 2 i = 3 i = 4
38
Richardson/Leapfrog
- Use two time step central differences
[( ) ] ( )
2 [( )] 2
O x x
u u u
x
u O t t
u u
t
u
n i
n i
n i
n
i
n i
n i
n
i
∂
−
∂
∂ (^) + −
- Result is explicit with second order
accuracy in time
( ) ( )
n i
n i
n i
n i
n i
n i
n i
n i
n i u u u u fu u u x
t u u 2 2 2
( )
2 1 1
Δ
αΔ = + + −
- However result is unstable for any f
and cannot be used
39
DuFort Frankel
- Modification of Richardson method to
provide stability
n
in second derivative by
average at time steps n+1 and n-
- Introduces another O[(Δt)
] error
⎥
⎥
⎦
⎤
⎢
⎢
⎣
⎡
Δ
Δ
−
( )
( ) ( ),( ), (^2) ( ) x
t O x t x
u T u u
t
u u
n i
n i
n i
n i
n i
n i
2 [( ) ]
1 1 2
u u u O t
n i
n i
n
i =^ + + Δ
[( ) ]
[( )]
O x
x
u u u
x
u
O t
t
u u
t
u
n i
n i
n i
n
i
n i
n i
n
i
40
DuFort Frankel
- Rearrange and introduce f = αΔt/(Δx)
( ) ( )
n i
n i
n i
n i
n i
n i
n i
n i
n i
n
i u u u u fu u u u
x
t
u u
( ) ( ) ( )
n
i
n
i
n
i
n
i
f u u f fu u
1 2 1 2 2
- Result is explicit for values at time n+
- Explicit start required to get first set of
values at time n-
Error in Explicit Solution of Conduction Equation
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
0.001 0.01 0.1 1 f = αΔ t/( Δ x)
2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Nx = 1000
Thermal Diffusivity = 1 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
Error in Crank-Nicholson Solution of Conduction
Equation
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
0.0001 0.001 0.01 0.1 1 10 100 1000 f = αΔ t/( Δ x)
2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Nx = 1000
Thermal Diffusivity = 1 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
Error in Fully Implicit Solution of Conduction Equation
0.
0.
0.
0.
0.
1
0.0001 0.001 0.01 0.1 1 10 100 1000 10000 100000 f = αΔ t/( Δ x)
2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Nx = 1000 Thermal Diffusivity = 1 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
Error in DuFort Frankel Solution of Conduction
Equation
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
0.0001 0.001 0.01 0.1 1 10 100 1000 f = αΔ t/( Δ x)
2
RMS Temperature Erro
Nx = 5 Nx = 10 Nx = 20 Nx = 50 Nx = 100 Nx = 200 Nx = 500 Nx = 1000
Thermal Diffusivity = 1 0 <= x <= (L = 1) T(x,0) = 1000 T(0,t) = T(L,t) = 0 End Time = 1
RMS Temperature Error by Execution Time
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E+
1.E+
1.E+
1.E+
1.E-02 1.E-01 1.E+00 1.E+01 1.E+02 1.E+03 1.E+04 1.E+ Execution Time (seconds)
RMS Temerature Error
Crank-Nicholson Explicit Fully Implicit DuFort-Frankel
46
This Week’s Homework
- Download assignment from web
- Find first and second derivative of sin x at x
= 1 for h = .1, .01, and .001 using second-
order central-difference expressions
- Repeat for x = 0.01 with h = .001 and.
- Do problems 5, 23, and 28 on pp 646-
of Hoffman (diffusion equation solutions)
- Computer assignment due after midterm
- Download and run program used to get charts
just shown for one of the suggested
assignments
47
Explicit Method Example
- How many values can you compute for
initial conditions below with f = 0.
( ) ( )
n
i
n
i
n
n i
i
n
i
n
i
n
i
u u u
u fu u fu +
??
??
??
? 70 80 90 100 90 80 70?
?
?
?
?
???
?????
80 90 95 90 80
88.75 (^) 92.588.
