1. Art and Diseño

Numerical Solutions to
Partial Differential Equations
Zhiping Li
LMAM and School of Mathematical Sciences
Peking University
A Model Problem in a 2D Box Region
Let us consider a model problem of parabolic equation:
ut = a (uxx + uyy ) ,
(x, y ) ∈ Ω, t > 0,
u(x, y , 0) = u 0 (x, y , 0),
(x, y ) ∈ Ω,
u(x, y , t) = 0,
(x, y ) ∈ ∂Ω, t > 0,
where a > 0 is a constant, Ω = (0, X ) × (0, Y ) ⊂ R2 .
1
For integers Nx ≥ 1 and Ny ≥ 1, let hx = 4x = XNx−1 and
hy = 4y = YNy−1 be the grid sizes in the x and y directions;
2
a uniform parallelopiped grid with the set of grid nodes
JΩ×R+ = {(xj , yk , tm ) : 0 ≤ j ≤ Nx , 0 ≤ k ≤ Ny , m ≥ 0},
where xj = j hx , yk = k hy , tm = mτ (τ > 0 time step size).
3
the space of grid functions
m
U = {Uj,k
= U(xj , yk , tm ) : 0 ≤ j ≤ Nx , 0 ≤ k ≤ Ny , m ≥ 0}.
Finite Difference Methods for Parabolic Equations
Explicit and Implicit Schemes in Box Regions
The Explicit Scheme in a 2D Box Region
The Forward Explicit Scheme in a 2D Box Region
The forward explicit scheme and its error equation
m+1
m
Uj,k
− Uj,k
m+1
ej,k
m
m
m
m
m
m
Uj+1,k
− 2Uj,k
+ Uj−1,k
Uj,k+1
− 2Uj,k
+ Uj,k−1
=a
+
.
τ
hx2
hy2
m
m
m
m
m
m
= [1 − 2(µx + µy )] ej,k
+µx ej+1,k
+ ej−1,k
+µy ej,k+1
+ ej,k−1
−Tj,k
τ,
where µx =
aτ
,
hx2
µy =
aτ
hy2
are the grid ratios in x and y directions.
The truncation error is
1
a
Tu(x, y , t) = utt (x, y , t)τ −
∂x4 u(x, y , t)hx2 + ∂y4 u(x, y , t)hy2
2
12
+ O(τ 2 + hx4 + hy4 ).
3 / 29
Finite Difference Methods for Parabolic Equations
Explicit and Implicit Schemes in Box Regions
The Explicit Scheme in a 2D Box Region
The Forward Explicit Scheme in a 2D Box Region
1
the condition for the maximum principle: µx + µy ≤ 12 ;
2
m = λm e i(αx xj +αy yk ) , α =
Fourier modes: Uj,k
x
α
3
h
amplification factor λα = 1 − 4 µx sin2
4
L2 stable if and only if µx + µy ≤ 1/2;
5
Convergence rate is O(τ + hx2 + hy2 ).
αx h x
2
lx π
X ,
αy =
+ µy sin2
ly π
Y ;
αy hy
2
i
;
4 / 29
The θ-Scheme in a 2D Box Region
"
#
"
#
δy2
δy2
δx2
δx2
m+1
m
= (1 − θ)a 2 + 2 Uj,k + θa 2 + 2 Uj,k
τ
hx
hy
hx
hy
m
m
m
m
m
m
Uj,k+1
− 2Uj,k
+ Uj,k−1
Uj+1,k − 2Uj,k
+ Uj−1,k
+
= (1 − θ)a
hx2
hy2
" m+1
#
m+1
m+1
m+1
m+1
m+1
Uj+1,k − 2Uj,k
+ Uj−1,k
Uj,k+1
− 2Uj,k
+ Uj,k−1
+ θa
+
.
hx2
hy2
m+1
m
Uj,k
− Uj,k
The error equation:
h i
m+1
m+1
m+1
m+1
m+1
[(1 + 2θ)(µx + µy )] ej,k
= θ µx ej+1,k
+ ej−1,k
+ µy ej,k+1
+ ej,k−1
+ (1 − θ) µx
m
+ [1 − 2(1 − θ)(µx + µy )] ej,k
m+ 1
m
m
m
m
ej+1,k
+ ej−1,k
+ µy ej,k+1
+ ej,k−1
− τ Tj,k 2 ,
Finite Difference Methods for Parabolic Equations
Explicit and Implicit Schemes in Box Regions
The θ-Scheme in a 2D Box Region
The θ-Scheme in a 2D Box Region
The truncation error
Tjm+∗
(
O(τ 2 + hx2 + hy2 ),
=
O(τ + hx2 + hy2 ),
if θ = 12 ,
if θ =
6 12 ,
1
the condition for the maximum principle:
2(µx + µy ) (1 − θ) ≤ 1;
2
m = λm e i(αx xj +αy yk ) , α =
Fourier modes: Uj,k
x
α
3
amplification factor
lx π
X ,
αy =
ly π
Y ;
h
i
α h
1 − 4(1 − θ) µx sin2 αx2hx + µy sin2 y2 y
h
i ;
λα =
α h
1 + 4θ µx sin2 αx2hx + µy sin2 y2 y
6 / 29
Finite Difference Methods for Parabolic Equations
Explicit and Implicit Schemes in Box Regions
The θ-Scheme in a 2D Box Region
The θ-Scheme in a 2D Box Region
4
for θ ≥ 1/2, unconditionally L2 stable;
5
for 0 ≤ θ < 1/2, L2 stable iff 2(1 − 2θ)(µx + µy ) ≤ 1;
6
the matrix of the linear system is still symmetric positive
definite and diagonal dominant when 2(µx + µy ) (1 − θ) ≤ 1,
however, each row has now up to 5 nonzero elements with a
band width of the order O(h−1 );
7
if solved by the Thompson method, the cost is O(h−1 ) times
of that of the explicit scheme;
8
in 3D, µx + µy ⇒ µx + µy + µz , and O(h−1 ) ⇒ O(h−2 ).
7 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
Alternative Approaches for Solving n-D Parabolic Equations
To reduce the computational cost, we may consider to apply highly
efficient iterative methods to solve the linear algebraic equations,
for example,
the preconditioned conjugate gradient method;
the multi-grid method;
etc..
Alternatively, to avoid the shortcoming of the implicit difference
schemes for high space dimensions, we may develop
the alternating direction implicit (ADI) schemes;
the locally one dimensional (LOD) schemes.
8 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
A 2D Alternating Direction Implicit (ADI) Scheme
A Fractional Steps 2D ADI Scheme by Peaceman and Rachford
1
in the odd fractional steps implicit in x and explicit in y , and
in even fractional steps implicit in y and explicit in x:
1
m+ 1
1 − µx δx2 Uj,k 2
2
1
m+1
1 − µy δy2 Uj,k
2
2
1
m
1 + µy δy2 Uj,k
,
2
1
m+ 1
=
1 + µx δx2 Uj,k 2 ,
2
=
numerical boundary conditions are easily imposed directly by
m+ 21
those of the original problem, since Uj,k
m+ 1
∼ uj,k 2 ;
9 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
A 2D Alternating Direction Implicit (ADI) Scheme
A Fractional Steps 2D ADI Scheme by Peaceman and Rachford
3
one step equivalent scheme
1
1
1
1
m+1
2
2
2
2
m
1 − µy δy Uj,k = 1 + µx δx
1 + µy δy Uj,k
.
1 − µx δx
2
2
2
2
4
Crank-Nicolson scheme
5
1
1
1
1
m+1
m
= 1 + µx δx2 + µy δy2 Uj,k
.
1 − µx δx2 − µy δy2 Uj,k
2
2
2
2
m+ 1
m+ 1
Since µx µy δx2 δy2 δt uj,k 2 = a2 τ 3 [uxxyyt ]j,k 2 + O(τ 5 + τ 3 (hx2 + hy2 )),
the truncation error of the scheme is O(τ 2 + hx2 + hy2 ).
10 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
A 2D Alternating Direction Implicit (ADI) Scheme
Stability of the 2D ADI Scheme by Peaceman and Rachford
1
m
i(αx xj +αy yk )
,
For the Fourier mode Uj,k
= λm
αe
λα =
αy hy
2
α h
2 µy sin2 y2 y
(1 − 2 µx sin2
)(1 − 2 µy sin2
)
(1 +
)(1 +
)
α x hx
2
2 µx sin2 αx2hx
,
the scheme is unconditionally L2 stable;
2
the two fractional steps can be equivalently written as
m+ 12
(1 + µx )Uj,k
µx
µy
m+ 1
m+ 21
m
m
+ Uj+1,k2 ,
Uj,k−1
+ Uj,k+1
Uj−1,k
+
2
2
1
1
µx
µy m+1
m+ 2
m+
m+1
Uj,k−1 + Uj,k+1
,
+
Uj−1,k
+ Uj+1,k2 +
2
2
m
= (1 − µy )Uj,k
+
m+ 21
m+1
(1 + µy )Uj,k
= (1 − µx )Uj,k
thus, the maximum principle holds if max{µx , µy } ≤ 1;
11 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
A 2D Alternating Direction Implicit (ADI) Scheme
Cost of the 2D ADI Scheme by Peaceman and Rachford
1
order the linear system by (Nx − 1) k + j and (Ny − 1) j + k in odd
and even steps, the corresponding matrixes are tridiagonal;
2
the computational cost: 3 times of that of the explicit scheme.
12 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
A 2D Alternating Direction Implicit (ADI) Scheme
The Idea of Peaceman and Rachford Doesn’t Work for 3D
In 3D, we can still construct fractional step scheme by
1
dividing each time step into 3 fractional steps;
2
introducing U m+ 3 and U m+ 3 at tm+ 1 and tm+ 2 ;
3
in the 3 fractional time steps, applying schemes which are in
turn implicit in x-, y - and z-direction and explicit in the other
2 directions respectively.
1
2
3
3
However, the equivalent one step scheme is, up to a higher order
term, the same as the θ-scheme with θ = 1/3, which has a local
truncation error O(τ + hx2 + hy2 ) instead of what we expect to have
O(τ 2 + hx2 + hy2 ) for an ADI method.
13 / 29
The Key Properties that Make an ADI Scheme Successful
1
Implicit only in one dimension in each fractional time step,
thus, if properly ordered, the matrix of the linear algebraic
equations is tridiagonal as well as diagonally dominant, hence
the computational cost is significantly reduced.
2
In the equivalent one step scheme, the implicit part is the
product of 1D implicit difference operators, the explicit part is
the product of 1D explicit difference operators (in a half time
step, + certain h.o.t. small perturbations), which guarantees
the scheme is unconditionally L2 stable;
3
The difference between the equivalent one step scheme and
the Crank-Nicolson scheme is a higher order term, which
guarantees that the local truncation error is O(τ 2 + h2 ).
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Alternate Direction Implicit (ADI) Schemes
Target Equivalent One Step Schemes of ADI Methods
For the 3D model problem, we aim to develop ADI finite difference
schemes which have a one time step equivalent scheme of the form
1
1
1
m+1
2
2
2
1 − µy δy
1 − µz δz Uj,k,l
1 − µx δx
2
2
2
1
1
1
2
2
2
m
= 1 + µx δx
1 + µy δy
1 + µz δz Uj,k,l
,
2
2
2
or 1
1
1
m+1
2
2
2
1 − µx δx
1 − µy δy
1 − µz δz Uj,k,l
2
2
2
1
1
1
2
2
2
m
= 1 + µx δx
1 + µy δy
1 + µz δz Uj,k,l
+ h.o.t.
2
2
2
15 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Alternate Direction Implicit (ADI) Schemes
An Extendable ADI Scheme Proposed by D’yakonov
For the 3D model problem, the ADI Scheme of D’yakonov is
µy 2 µz m
µx 2 µx 2 m+∗ δx Uj,k,l = 1 +
δx 1 +
δy 1 + δz2 Uj,k,l
,
1−
2 2
2
2
µy 2
m+∗∗
m+∗
1−
δy Uj,k,l
= Uj,k,l
,
2
µz
m+1
m+∗∗
1 − δz2 Uj,k,l
= Uj,k,l
.
2
The scheme obviously has all the key properties. However, there is
one thing need to be taken care of: U m+∗ , U m+∗∗ are not
approximate solutions, their boundary conditions may not be
directly derived from that of the original problem.
16 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Alternate Direction Implicit (ADI) Schemes
Numerical Boundary Conditions for U m+∗ , U m+∗∗
m+1
m+1
m+1
m+1
Since Uj,k,l
∼ uj,k,l
, Uj,k,l
= uj,k,l
on the boundary (j = 0, Nx , or
k = 0, Ny , or l = 0, Nz ). Thus, the boundary condition
m+∗∗
Uj,k,l
: j = 0, Nx , 0 ≤ k ≤ Ny , 0 < l < Nz , or k = 0, Ny ; 0 ≤ j ≤ Nx , 0 < l < Nz
m+∗∗
can be obtained from Uj,k,l
= 1−
µz 2 m+1
2 δz Uj,k,l .
The boundary
n condition
o
m+∗
Uj,k,l
: j = 0, Nx , 0 < k < Ny , 0 < l < Nz .
m+∗∗
µ
m+∗
can be obtained from Uj,k,l
= 1 − 2y δy2 Uj,k,l
and the
m+∗∗
boundary values for Uj,k,l
on the nodes j = 0, Nx .
17 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Alternate Direction Implicit (ADI) Schemes
An Extendable ADI Scheme Proposed by Douglas and Rachford
For the 3D model problem, the Scheme is of the form
1
m+1∗
m+1∗
m
m
m
m
+ Uj,k,l
Uj,k,l
= Uj,k,l
+ µx δx2 Uj,k,l
+ µy δy2 Uj,k,l
+ µz δz2 Uj,k,l
,
2
1
m+1∗∗
m+1∗∗
m+1∗
m
m
+ Uj,k,l
Uj,k,l
= Uj,k,l
+ µy δy2 Uj,k,l
− µy δy2 Uj,k,l
,
2
1
m+1
m+1
m+1∗∗
m
m
+Uj,k,l
− µz δz2 Uj,k,l
.
Uj,k,l
= Uj,k,l
+ µz δz2 Uj,k,l
2
Since U m+1∗ and U m+1∗∗ are approximate solutions at tm+1 , their
boundary conditions can be directly derived from that of the
original problem.The scheme has all the key properties with the
m .
h.o.t.= − 14 µx µy µz δx2 δy2 δz2 Uj,k,l
18 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Alternate Direction Implicit (ADI) Schemes
The Idea of the ADI Scheme of Douglas and Rachford
The construction of the scheme can be viewed as a predictcorrection process.
In the 1st fractional step, Crank-Nicolson to x, explicit to y , z.
In the 2nd fractional step, to improve stability and accuracy in y ,
the y -direction is corrected by the Crank-Nicolson scheme:
1
m+1∗
m+1∗∗
m
m
+ Uj,k,l
Uj,k,l
= Uj,k,l
+ µx δx2 Uj,k,l
2
1
m+1∗∗
m
m
+ µy δy2 Uj,k,l
+ Uj,k,l
+ µz δz2 Uj,k,l
.
2
Note, this plus the 1st equation gives the second in the scheme.
19 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Alternate Direction Implicit (ADI) Schemes
The Idea of the ADI Scheme of Douglas and Rachford
In the 3rd, z-direction is corrected by the Crank-Nicolson scheme:
1
m+1∗
m+1
m
m
+ Uj,k,l
Uj,k,l
= Uj,k,l
+ µx δx2 Uj,k,l
2
1
1
m+1∗∗
m+1
m
m
2
+ Uj,k,l
+ µz δz2 Uj,k,l
+ Uj,k,l
.
+ µy δy Uj,k,l
2
2
Note, the 3rd equation in the scheme is simply the difference of
the above two equations.
20 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Locally One Dimensional (LOD) Schemes
An Extendable Locally One Dimensional (LOD) Scheme
1
In the ith fractional step, the problem is treated as a one
dimensional problem of the ith dimension.
2
1D Crank-Nicolson scheme is applied to each fractional step.
For the 3D
the LOD
is given as
model problem,
scheme 3
1
m+∗
1 − µx δx2 Uj,k,l
2
1
m+∗∗
1 − µy δy2 Uj,k,l
2
1
m+1
1 − µz δz2 Uj,k,l
2
4
=
=
=
1+
1+
1+
1
m
µx δx2 Uj,k,l
,
2
1
m+∗
µy δy2 Uj,k,l
,
2
1
m+∗∗
.
µz δz2 Uj,k,l
2
The scheme has all of the key properties.
21 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Locally One Dimensional (LOD) Schemes
Impose Boundary Conditions for the LOD Scheme
U m+∗ , U m+∗∗ are nonphysical, their boundary conditions must be
handled with special care.
1
multiplying (1 − 12 µz δz2 ) on the 3rd equation leads to
1 2 4
1 2 4
m+∗∗
m+1
2
,
1 − µz δz Uj,k,l = 1 − µz δz + µz δz Uj,k,l
4
4
Hence, by omitting the O(τ 2 ) order terms, the boundary
condition for U m+∗∗ can be given by:
m+1
m+∗∗
Uj,k,l
= 1 − µz δz2 Uj,k,l
,
k = 0, Ny , 0 ≤ j ≤ Nx , 0 < l < Nz ,
j = 0, Nx , 0 ≤ k ≤ Ny , 0 < l < Nz .
22 / 29
Finite Difference Methods for Parabolic Equations
2D and 3D ADI and LOD Schemes
3D Locally One Dimensional (LOD) Schemes
Impose Boundary Conditions for the LOD Scheme
2
Similarly, multiplying (1 − 12 µy δy2 ) on the 2nd equation,
omitting the O(τ 2 ) and higher order terms, the boundary
condition for U m+∗ can be given by:
m+∗∗
m+∗
Uj,k,l
= 1 − µy δy2 Uj,k,l
,
j = 0, Nx , 0 < k < Ny , 0 < l < Nz .
23 / 29
Finite Difference Methods for Parabolic Equations
FDM for General Parabolic Problems in High Dimensions
General Domain and Boundary Conditions
General Domain and Boundary Conditions
1
Approximate boundary conditions can be established with the
similar methods used in § 1.3.4 for elliptic problems.
2
Approximate boundary conditions can further restrict the
stability condition for explicit schemes, thus implicit schemes
are even more preferable.
3
ADI and LOD schemes are in more complicated forms.
4
Some special regions with non-planar boundaries can be
transformed into box regions by introducing curved coordinate
systems, say, the polar coordinate system can be used for a
circular region, the cylindrical coordinate system can be used
for cylindrical regions, etc...
24 / 29
Finite Difference Methods for Parabolic Equations
FDM for General Parabolic Problems in High Dimensions
Variable-coefficient Equations
Variable-coefficient and nonlinear Equations
1
ADI and LOD schemes can also be extended to high
dimensional variable-coefficient linear, and even certain
nonlinear problems.
2
In such cases, the difference operators (1 ± 12 µx δx2 ) and
(1 ± 12 µy δy2 ) are generally noncommutative, i.e.
1
1
1
1
2
2
2
2
1 ± µx δx
1 ± µy δy 6= 1 ± µy δy
1 ± µx δ x ,
2
2
2
2
and this will introduce the so called splitting error.
3
Boundary conditions are more difficult to handle.
25 / 29
Finite Difference Methods for Parabolic Equations
FDM for General Parabolic Problems in High Dimensions
Variable-coefficient Equations
Implicit, Semi-implicit Schemes
For nonlinear problems, nonlinear algebraic equations derived from
implicit schemes need to be solved with iterative methods, such as
semi-implicit schemes.
The idea is to
apply an implicit scheme only to the principal linear part of
the nonlinear equation;
approximate the residual nonlinear part by an explicit scheme.
26 / 29
Finite Difference Methods for Parabolic Equations
FDM for General Parabolic Problems in High Dimensions
Variable-coefficient Equations
Implicit, Semi-implicit Schemes and Elliptic Solver
To apply the implicit schemes to linear parabolic problems and the
semi-implicit schemes to nonlinear parabolic problems, it is a
sequence of linear algebraic systems corresponding to
certain elliptic problems we actually need to solve in
the end.
So, Fast solvers for elliptic problems play important roles in solving
parabolic problems.
27 / 29
Finite Difference Methods for Parabolic Equations
FDM for General Parabolic Problems in High Dimensions
Variable-coefficient Equations
Asymptotic Analysis and Extrapolation Methods
The asymptotic analysis and extrapolation methods introduced in
§ 1.5 for elliptic problems can also be extended to analyze the
finite difference approximation error for parabolic problems.
In particular, the error bounds can be estimated by the numerical
solutions obtained on grids with different mesh sizes.
For example, for the explicit scheme of the heat equation, let
m = u m + O(h2 ) + O(h4 ),
τ = µh2 , µ ≤ 1/2, then, we have U(h)j
(h)j
(1)m
Uj
,
4m
m
4U(h/2)2j
− U(h)j
3
m
= u(h)j
+ O(h4 ).
28 / 29
Thank You!
SK 2µ20, 21