...

AEM2-Ch07.pdf

by user

on
Category: Documents
2

views

Report

Comments

Transcript

AEM2-Ch07.pdf
Chapter 7
Numerical Methods for Partial Differential
Equations
1
Introduction
Many problems encountered in applications caanot be solved using any of the analytical
methods previously described. It is then necessary to resort to approximate solution methods. The most commonly used approximation methods for the solution of partial differential
equations are the finite difference method, the finite volume method and the finite element
method. In these methods one represents the solution of the problem in a finite dimensional
space of points called nodes. Appropriate discrete forms of the governing equations are used
to obtain computational algorithms that allow calculation of the approximate values of the
solution at the nodal locations.
Many problems in engineering and science are formulated as partial differential equations.
As previously noted, the following types of partial differential equations are of frequent
occurence in applications.
The elliptic equation known as the Laplace’s equation
∂ 2u ∂ 2u
+
=0
∂x2 ∂y 2
satisfied by the function u(x, y), describes all sorts of steady state problems. Laplace’s
equation arises in the study of steady state temperature distribution in an isotropic medium,
the analysis of steady state current distribution in homogeneous conductors and in the
analysis of the velocity potential due to an irrotational incompressible fluid flow.
Poisson’s equation
∂ 2u ∂ 2u
+
= f(x, y)
∂x2 ∂y 2
satisfied by the function u(x, y), represents the problem of conduction of heat inside a planar
domain undergoing internal heat generation. Poisson equation is encountered in the theory
1
of torsion, the theory of slow-moving incompressible viscous flow and the inverse-square
law theories of electricity, magnetism and gravitation in the presence of charge densities,
magnetic poles and mass densities.
The parabolic equation known as the heat or diffusion equation
∂u
∂ 2u
=α 2
∂t
∂x
satisfied by the function u(x, t), describes for example the transport of energy or diffusant
through a conducting or diffusing medium.
Finally the hyperbolic equation called the wave equation
2
∂ 2u
2∂ u
=
c
∂t2
∂x2
satisfied by the function u(x, t), describing the instantanoeus shape of a vibrating string.
In each case, the above equations must be solved for u(x, y) or u(x, t) subject to specified
boundary and initial conditions.
2
2.1
Numerical Methods for Elliptic Problems
Finite Difference Method
Consider the case of Poisson’s equation for the unknown function u(x, y),
∂ 2u ∂ 2u
+
= f(x, y)
∂x2 ∂y 2
on (x, y) in the interior R = {(x, y) : a < x < b, c < y < d} of a planar region and subject to
u(x, y) = g(x, y)
along the boundary S of the region. As long as f and g are continuous, one can show that
a unique solution exists.
A simple approach to the numerical solution of the above problem consists in partitioning
the intervals [a, b] and [c, d] respectively into n and m subintervals with step sizes ∆x = h =
(b−a)/n and ∆y = k = (d−c)/m so that the node or mesh point (xi, yj ) is at xi = a+(i−1)h
for i = 1, .., N and yj = c + (j − 1)k for j = 1, .., M. Here, N = n + 1 and M = m + 1.
In practice, a distinction must be made between interior nodes (i.e. (xi , yj ) for i =
2, ..., N − 1 and j = 2, ..., M − 1) and the boundary nodes (represented by nodes with i = 1,
i = N, j = 1, or j = M). In the finite difference method, all interior nodes are handled by the
same equation. Special equations are typically required for the boundary nodes depending
of the boundary conditions applied in each case.
2
Using now centered finite difference approximations for the partial derivatives, the following second order accurate system of simultaneous algebraic equations is obtained for all
interior nodes
h
h
2[( )2 + 1]wi,j − (wi+1,j − wi−1,j ) − (wi,j+1 − wi,j−1 )( )2 = −h2f(x, y)
k
k
This expression is often called the five point formula.
If the non-homogeneous term is f(x, y) = 0 (Laplace’s equation) and if the mesh spacings
are chosen to be the same (i.e. h = k) the following particularly simple system of equations
is obtained for the interior nodes:
1
wi,j = [wi+1,j + wi−1,j + wi,j+1 + wi,j−1 ]
4
The complete set of equations for all nodes in the mesh constitutes a system of linear
algebraic equations. The system may be solved by a direct method (Gauss elimination) or
more commonly, by an iterative method. In an iterative method one assumes an initial guess
old
new
for the solution wi,j
and computes an improved guess wi,j
using the discretization formula.
The process is then repeated until the computed values at all grid points change little with
every subsequent iteration (convergence of the iterations). Typically, at every iteration once
searches the entire mesh for the largest relative difference between iterates, i.e.
ǫM AX = MAX(ǫi,j ) = MAX(
new
old
|wi,j
− wi,j
|
|)
old
|wref
old
where wref
is a suitably selected reference value. Iterations stop once a predetermined level
of tolerance is reached, i.e.
ǫM AX ≤ ǫT OL
In the Jacobi iteration procedure, the iteration formula is
1 old
new
old
old
old
wi,j
= [wi+1,j
+ wi−1,j
+ wi,j+1
+ wi,j−1
]
4
In the Gauss-Seidel iteration procedure one takes advantage of the fact that nodes are
always visited in sequence (say, from laft to right and bottom to top) and uses the latest
values stored in the memory of the computer, as soon as these become available. The
iteration formula in this case is
1 old
new
new
old
new
wi,j
= [wi+1,j
+ wi−1,j
+ wi,j+1
+ wi,j−1
]
4
In the Successive Over-relaxation (SOR) iteration procedure one improves on the G-S
procedure by relaxing the value of the iterated variable. The iteration formula in this case is
ω old
new
old
new
old
new
old
wi,j
= wi,j
+ [wi+1,j
+ wi−1,j
+ wi,j+1
+ wi,j−1
− 4wi,j
]
4
3
where ω (1 < ω < 2) is the relaxation factor. Note that the SOR formula reduces to the
Gauss-Seidel scheme when ω = 1.
Exercise. Consider steady state heat conduction in a thin square plate of edge = 100
cm. The edges x = 0, x = l and y = 0 are maintained at u = 0 while at the edge y = d
u(x, d) = 100. No heat flow along the z direction perpendicular to the plate. The required
temperature u(x, y) satisfies
∂ 2u ∂ 2u
+
=0
∂x2 ∂y 2
This problem is readily solved in closed form by separation of variables yielding
nπ
y)
nπ sinh( 100
u(x, y) =
cn sin(
x)
100 sinh(nπ)
n=1
∞
X
with
cn =
(
0
400
nπ
n even
n odd
Solve the above problem numerically using the finite difference method and compare your
results against the analytical solution.
The five point formulae above are also applicable at the boundary nodes but they must be
supplemented by specially formulated discrete equations representing the boundary conditions of the problem in order to obtain the specific finite difference formulae for the boundary
nodes.
If Dirichlet boundary conditions are imposed on the boundary nodes then one has
w1,j = g(x1, yj ) j = 1, .., M
wN,j = g(xN , yj ) j = 1, .., M
wi,1 = g(xi , y1) i = 1, 2, .., N
wi,M = g(xi , yM ) i = 1, 2, .., N
To obtain a banded matrix the mesh points must be relabeled sequentially from left to
right and from top to bottom. The resulting system can be solved by Gaussian elimination
if n and m are small and by iterative methods (e.g. SOR iteration) when they are large.
One is clearly interested in determining the quality of the approximation computed by the
finite difference method. This is measured by the error. By using Taylor series expansions,
the following error bound for the numerical solution is obtained for the case when ∆x = h =
∆y = k = h
h2
|wi,j − u(xi, yj )| ≤ (Mxxxx + Myyyy )
96
where Mxxxx and Myyyy are bounds for uxxxx and uyyyy .
4
2.2
Finite Volume Method
An alternative discretization method is based on the idea of regarding the computation
domain as subdivided into a collection of finite volumes. In this view, each finite volume is
represented by an area in 2D and a volume in 3D. A node, located inside each finite volume
is the locus of computational values. In rectangular cartesian coordinates in 2D the simplest
finite volumes are rectangles. For each node, the rectangle faces are formed by drawing
perpendiculars through the midpoints between contiguous nodes. Discretization equations
are obtained by integrating the original partial differential equation over the span of each
finite volume. The method is easily extended to nonlinear problems.
Consider for instance the following problem consisting of determining u(x, y) in x ∈ [a, b],
y ∈ [c, d] such that
∂ ∂u
∂ ∂u
(a ) +
(a ) = f
∂x ∂x
∂y ∂y
subject to specified conditions at the boundary. The problem represents steady state heat
conduction in a solid with position dependent conductivity a where energy is being internally
generated at the rate f, per unit volume.
The typical control volume has dimensions ∆x × ∆y × 1 and the discrete equation is
obtained by perfoming an energy balance on the finite volume. The net energy input is
obtained by integrating the fluxes a∂u/∂x and a∂u/∂y into and out from the finite volume
in the x and y directions, respectively. Since energy is conserved, the result must equal the
rate of energy generation inside the volume which is given by the volume integral of f.
A notation from carthography is used to simplify the discrete formulation. In this notation, the node representing the finite volume is called P and the neighboring nodes along the
four coordinate directions are called N, S, E and W , respectively. The following discrete
equation is obtained for node P ,
ae (wE − wP ) + aw (wP − wW ) an (wN − wP ) + as (wP − wS )
+
= fP
∆x2
∆y 2
where the subscripts e, w, n, s are used to indicate that the values of a are to be evaluated at
the corresponding finite volume faces. Note that if ae = aw = an = as = a = 1, ∆x = ∆y = h
and fP = 0 the above reduces to the simple form
1
wP = [wE + wW + wN + wS ]
4
which is identical to the expression obtained before using the method of finite differences.
As in the case of finite differences, care is required to properly deal with the boundary
conditions at the boundary nodes, specially when derivative conditions are involved.
When the boundaries of the computational domain do not coincide with coordinate lines
additional care must be used for the implementation of boundary conditions. Finite difference
5
or finite volume formulae can be derived by ad-hoc methods, by conventional formulae if the
irregular boundaries are represented by staicases approximating the actual boundaries or by
coordinate transformation methods.
2.3
Iterative Solution of the Algebraic System
The system of algebraic equations is readily solved using iterative methods. Specifically, in
the Jacobi method for the numerical solution of Laplace’s equation in a uniform mesh one
assumes aa solution for the interior nodes and then computes an improved approximation
using the five point formula, i.e.
1 o
o
o
o
+ wi−1,j
+ wi,j+1
+ wi,j−1
]
wi,j = [wi+1,j
4
0
where wi,j
are the initial guesses values and wi,j is the improved guess generated by iteration.
Since nodes [i − 1, j] and [i, j − 1] would have been already visited when computing the
improved guess for node [i, j], it is wise to use the improved values as soon as they become
available in the memory of the computer. This is Gauss-Seidel method and is expressed as
1 o
o
wi,j = [wi+1,j
+ wi−1,j + wi,j+1
+ wi,j−1 ]
4
Finally, if the size of the corrections introduced by iteration is artifically enhanced by
means of an overrelaxation factor one obtains the successive overrelaxation (SOR) method,
i.e.
ω o
o
o
o
+ wi−1,j + wi,j+1
+ wi,j−1 + 4wi,j
]
wi,j = wi,j
+ [wi+1,j
4
where ω > 1 is the overrelaxation factor. Note that this formula reduces to G-S when ω = 1.
The iteration process stops whenever the maximum difference found between two contiguous iterates falls below a specified limit.
2.4
Finite Element Method
The finite element method is an alternative discrete representation obtained from the variational formulation of the original problem. Specifically, for the problem
∂ 2u ∂ 2u
+ 2 = ∇2 u = f
2
∂x
∂y
inside the domain Ω = {x|x ∈ [a, b]} and {y|y ∈ [c, d]}, and subject to u = 0 at the boundary,
the associated functional I is given by
I(v) =
1
[ |∇v|2 + fv]dxdy
Ω 2
Z
6
The desired solution is the function u ∈ v which minimizes the functional, i.e.
δI(u) = 0
and this function comes from the set of functions which satisfy the stated boundary conditions.
Discretization in the finite element method is also obtained by subdiving the computational domain into subdomains (finite elements), commonly triangles, rectangles, tetrahedra
or rectangular parallelepipeds (bricks). The vertices in each of these geometries are called
nodes. The methods seeks an approximation inside each element of the form
w(x, y) =
m
X
ci φ i
i=1
where the φi are linearly independent piecewise polynomials and the ci are constants with
∂I
= 0). The constants ci actually
respect to which the functional is minimized (i.e. ∂c
i
represent the values of u at the nodal locations.
Minimization of I with respect to the ci ’s produces a set of simultaneous linear algebraic
equations of the form
Ac = b
Solution of the above system using standard techniques yields the desired approximation.
Consider as a second example the following partial differential equation in a plane region
D
∂
∂u
∂
∂u
(p(x, y) ) +
(q(x, y) ) + r(x, y)u(x, y) = f(x, y)
∂x
∂x
∂y
∂y
subject to the boundary conditions
u(x, y) = g(x, y)
along portion S1 of the boundary S of D and along the remaining portion S2 of the boundary
p
∂u
∂u
cos θ + q
cos(π/2 − θ) + g1 (x, y)u = g2(x, y)
∂x
∂y
where θ is the angle of the outward normal.
Now, in D let p and q be C 1, r, f, g1 and g2 be all C0, let also p > 0, q > 0 and r ≤ 0.
Then the solution of the PDE above is the twice differentiable function w(x, y) satisfying
w = g which minimizes the functional
I[w] =
1
∂w
∂w
{ [p(x, y)( )2 + q(x, y)( )2 − r(x, y)w2 ] + f(x, y)w}dxdy +
∂x
∂y
D 2
Z
1
{−g2 w + g1 w2 }dS
2
S2
Z Z
7
The search for w is performed in the FEM from among the (smaller) set of functions which
are linear combinations of appropriately chosen basis function just as in the Rayleigh-RitzGalerkin method.
As before, the region D is first subdivided into a finite number of regions with simple
shapes (e.g. triangles or rectangles). The vertices of the resulting network are called nodes.
Piecewise polynomials basis functions for triangle shaped elements are
φi (x, y) = a + bx + cy
The method then seeks to determine constants γi such that
w(x, y) ≈ φ =
minimizes I[w] ≈ I[
Pm
i=1
m
X
γi φi (x, y)
i=1
γi φi (x, y)]. The condition for a minimum is
∂I
=0
∂γj
for each j = 1, 2, ..., n. Note that constants γn+1 , γn+2 , ..., γm are determined by the boundary
condition φ = g.
Applying the condition for a minimum and taking advantage of the simple form of φi
leads to the system
Ac = b
where the entries of A and b involve double integrals over D and path integrals over S2. The
resulting system can then be solved using standard methods.
Exercise. A cantilever beam of length L has rectangular cross section with height h = 2a
and width w = 2b. A point load of magnitude P is applied on the free end of the beam
and perpendicular to its length. Recall that the moment of inertia of the beam cross section
with respect to a centroidal axis is given by
I=
1
wh3
12
Introducing rectangular cartesian coordinates with the x−y plane parallel to the cross section
of the beam and the z axis along its length, and in terms of the stress function φ(x, y) given
by
τxz =
∂φ
∂y
and
τyz = −
∂φ
∂x
8
the problem becomes
∂ 2φ ∂ 2φ
ν Py
+ 2 =
2
∂x
∂y
1+ν I
which is Poisson’s equation and it also represents the deflection of a uniformly stretched
rectangular membrane under a continuous load.
The above problem can be solved in closed form to give
φ(x, y) = −
X n=∞
X (−1)m+n−1 cos[(2m + 1)πx/2a] sin(nπy/b)
ν P 8b3 m=∞
1 + ν I π 4 m=0 n=1
(2m + 1)n[(2m + 1)2 (b2/4a2 ) + n2 ]
Solve the above problem numerically using the finite element method and compare your
results against the analytical solution.
Exercise. A thin rectangular plate of edges a and b and thickness h is constrained along
all four edges and is loaded by hydrostatic pressure q acting uniformly on one side. Compute
the elastic deflection of a plate of sides a = 1 m, b = 0.75 m, and thickness h = 0.025
m, uniformly loaded by pressure q = 5 × 105 (N/m2 ), using the finite element method and
compare against the analytical solution given by
w(x, y) =
∞
(−1)(m−1)/2
mπx
αm tanh(αm ) + 2
mπy
4qa4 X
cos(
)[1 −
cosh(
)
5
5
π D m=1,3,5,...
m
a
2 cosh(αm )
a
+
1
mπy
mπy
sinh(
)]
2 cosh(αm ) a
a
where
αm =
mπb
2a
and
D=
Eh3
12(1 − ν 2 )
is the flexural rigidity of the plate. Assume E = 2 × 1011 (N/m2 ), ν = 0.33.
3
Parabolic Problems: Finite Difference Methods
The finite difference method is readily applied to the solution of parabolic problems. In this
case a mesh of nodes along the time coordinate is introduced in addition to the spatial mesh.
The various derivatives in the parabolic equation are then approximated by finite difference
formulae. The resulting system of algebraic equations is then solved at each time level in
order to obtain the approximate values of the solution at the spatial nodes. Depending
9
on how one approaches the problem along the time coordinate, three different schemes are
obtained: the explicit scheme, the implicit scheme and the semi-inplicit scheme. These are
described in detail below.
As described before, the simplest parabolic problem consists of determining the function
u(x, t) satisfying the parabolic equation
∂u
∂ 2u
=α 2
∂t
∂x
subject to suitable initial and boundary conditions.
3.1
Explicit Scheme
As mentioned above, in the finite difference method the space domain is subdivided into
equal size subintervals to give the set of mesh points (in space) xi = (i − 1)∆x = (i − 1)h
with i = 1, 2, ..., N while a mesh of points in time tj = (j−1)∆t = (j−1)k with j = 1, 2, ..., M
is also created.
If a forward difference approximation is now used for the time derivative and a centered
difference formula is used for the spatial derivative the following first order accurate in time,
second order accurate in space, explicit scheme is obtained
wi+1,j − 2wi,j + wi−1,j
wi,j+1 − wi,j
=α
k
h2
As the name implies, this formulae can be solved explicitly to give the values of w for all
points on the space mesh at a subsequent time level given the values of w at all mesh points
at the immediately preceeding time level, i.e.
wi,j+1 = (1 −
2αk
k
)wi,j + α 2 (wi+1,j + wi−1,j ) = (1 − 2αν)wi,j + αν(wi+1,j + wi−1,j )
2
h
h
where ν = k/h2 . The above formula can also be written using matrix notation as w(j) =
Aw(j−1).
Note that each equation in the above set uses the values of w at three contiguous nodes
(i − 1, i, i + 1) for time level j to compute the value of wi,j+1 .
The explicit method, however, is only conditionally stable. This means that only for
certain predetermined values of the mesh spacings h and k is the truncation error damped
and thus unable to increase as the calculation progresses. For stability, it is required that the
spectral radius ρ(A) ≤ 1. As a result, for stability of the explicit scheme the mesh spacing
must satisfy the Courant-Friedrichs-Lewy (CFL) restriction
k
1
=ν≤
2
h
2α
10
Stability is frequently investigated by expanding the numerical solution in Fourier modes,
i.e.
wi,j = (λ)j exp[ıkxi∆x]
√
where ı = −1, λ is the error amplification factor and kx is the error wave number, and
substituting into the finite difference formula to obtain a relationship for the amplification
factor which in turn yields the CFL restriction.
The truncation error of the scheme, T (x, t) is given by
T (x, y) =
u(xi, tj+1 ) − u(xi, tj )
u(xi+1 , tj ) − 2u(xi, tj ) + u(xi−1 , tj )
−α
=
k
h2
1
1
= utt(x, η)k − uxxxx(ξ, t)h2
2
12
where x ≤ ξ ≤ x + h and t ≤ η ≤ t + k. An upper bound for the truncation error is
1
1
T (x, y) ≤ k[Mtt + Mxxxx ]
2
6ν
where the Ms are bounds for the corresponding derivatives. The explicit scheme is consistentwith the original PDE since
T →0
as k → 0.
The explicit finite difference scheme is called convergent if the computed solution at any
point approaches the exact value as the mesh is refined, i.e. if
wi,j → u(xi , tj )
1
as h → 0, k → 0. It can be shown that the explicit scheme is convergent as long as ν ≤ 2α
.
Exercise As mentioned before, the transport of diffusing species and of energy through
and out of a solid can be described by a parabolic partial differential equation equation (often
called the heat equation or the diffusion equation). Consider the following specific example
consisting of determining u(x, t) such that
∂u
∂ 2u
(x, t) = α 2 (x, t)
∂t
∂x
on 0 < x < l, 0 < t and subject to the conditions
u(0, t) = u(l, t) = 0,
t>0
u(x, 0) = f(x), 0 ≤ x ≤ l
11
Here u(x, t) is the concentration of diffusing substance (in mass transfer problems) or the
temperature (in energy transfer problems), at position x and time t, f(x) is the initial
distribution and α is the diffusivity. This problem is readily solved using separation of
variables. Assuming α = 1, the solution is
u(x, t) =
∞
X
am sin(mπx) exp[−(mπ)2t]
m=1
where
am = 2
Z
0
1
f(x) sin(mπx)dx
Solve the above problem using the method of finite differences with an explicit scheme.
Assume α = 1 and f(x) = 2x for 0 ≤ x ≤ 21 and f(x) = 2 − 2x for 12 ≤ x ≤ 1. Use m = 20
and time steps just above and just below that given by the CFL condition. Compare against
the analytical solution.
3.2
Implicit Scheme
An unconditionally stable scheme results for linear problems if backward differences in time
are used instead, i.e.
wi,j − wi,j−1
wi+1,j − 2wi,j + wi−1,j
=α
k
h2
which is equivalent to
wi,j+1 − wi,j
wi+1,j+1 − 2wi,j+1 + wi−1,j+1
=α
k
h2
This practice produces after slight rearrangement the implicit method given by the formula
−
αk
αk
αk
wi−1,j+1 + [1 + 2( 2 )]wi,j+1 − 2 wi+1,j+1 = wi,j
2
h
h
h
For a given time level j + 1, the above can be represented as
−ai wi−1,j+1 + bi wi,j+1 − ci wi+1,j+1 = di
for each of the xi nodes in the spatial mesh. Note that each equation in the above set
involves the values of w at three contiguous nodes (i − 1, i, i + 1) for time level j + 1 and wi,j
to compute the value of wi,j+1 , i.e. this is a system of simultaneous, interlinked algebraic
equations. The system of equations is clearly associated with a tridiagonal matrix and is
written in matrix notation as Aw(j) = w(j−1) .
12
The implicit scheme is unconditionally stable as can be shown by representing the numerical solution in terms of its Fourier modes. However, although the scheme is consistent, the
trunctation error is only O(k). In any case, since consistency and stability imply convergence
(Lax equivalence theorem), the implicit scheme is convergent.
Implicit Finite Difference Heat Equation Algorithm.
• Give initial condition f(x, y), boundary conditions, number of subintervals in x and y
and time step ∆t.
• Create coefficients for the tridiagonal matrix.
• Solve tridiagonal matrix.
• Move to the next time step.
Exercise. Same as exercise above but using the implicit scheme.
3.3
The Semi-Implicit (Crank-Nicolson) Scheme
The explicit and implicit schemes can be combined to produce a semi-implicit scheme with
improved features. The schemes can be combined by simple weighing
wi,j+1 − wi,j α
wi+1,j − 2wi,j + wi−1,j
wi+1,j+1 − 2wi,j+1 + wi−1,j+1
− [(1 − θ)
+
θ
]=0
k
2
h2
h2
where 0 ≤ θ ≤ 1 is the weighing factor. Note that θ = 0 gives the explicit scheme and θ = 1
yields the implicit scheme.
A stability analysis of this generalized scheme can be performed using Fourier modes and
the result is that an unconditionally stable scheme is obtained as long as 21 ≤ θ ≤ 1.
A second order accurate in time and space scheme can be obtained by averaging the
forward and backward differences approximations (i.e. θ = 12 to produce the Crank-Nicolson
method
wi,j+1 − wi,j α wi+1,j − 2wi,j + wi−1,j wi+1,j+1 − 2wi,j+1 + wi−1,j+1
− [
+
]=0
k
2
h2
h2
where the associated matrix A is positive definite, diagonally dominant and tridiagonal.
Either the Crout scheme or SOR iteration can be used to solve this system. Note that each
equation in the above set uses the values of w at three contiguous nodes (i − 1, i, i + 1) for
time level j and the same three nodes at time level j + 1 to compute the value of wi,j+1 .
Crank-Nicolson Finite Difference Algorithm.
• Give initial condition f(x, y), boundary conditions, number of subintervals in x and y
and time step ∆t.
• Create coefficients for the tridiagonal matrix.
13
• Solve tridiagonal matrix.
• Move to the next time step.
Exercise. Same as exercise above but using the C-N scheme.
3.4
More General Boundary Conditions
Special care is require to handle more general boundary conditions. As an example consider
the situation where the following condition applies at x = 0
∂u
= βu + g
∂x
where β and g can be functions of time. To derive a discrete equation for the node located
at x = 0 one introduces a fictitious node located at x = −h and uses centered differences to
write
w+h,j − w−h,j
w2,j − w−2,j
=
= βw1,j + g
2h
2h
The fictitious node value w−h,j is then eliminated by combining the above with the finite
difference formula and scheme being used for the interior nodes (explicit, implicit or semiimplicit) which is assumed valid up to x = −h. The result is an algebraic equation for
w1,j .
3.5
One-dimensional Problems in Polar Coordinates and Boundary Conditions
When a multidimensional problem exhibits cylindrical or spherical symmetry it can be represented as a one dimensional problem in polar coordinates, i.e. for a cylindrical coordinate
system
∂u
1 ∂ ∂u
∂ 2u 1 ∂u
1
1
=
(r ) = 2 +
= ut = (rur )r = urr + ur
∂t
r ∂r ∂r
∂r
r ∂r
r
r
and for a spherical coordinate system
∂u
1 ∂
∂u
∂ 2u 2 ∂u
1
2
= 2 (r2 ) = 2 +
= ut = 2 (r2 ur )r = urr + ur
∂t
r ∂r
∂r
∂r
r ∂r
r
r
The above two equations can be represented by the single formula
γ
ut = urr + ur
r
where γ = 1 for cylindrical and γ = 2 for spherical coordinates.
14
The corresponding second order accurate finite difference formula using an explicit scheme
and a uniform mesh in space are then readily obtained. Specifically, in cylindrical coordinates
1
wi,j+1 − wi,j
=
(r 1 wi+1,j − 2ri wi,j + ri− 1 wi−1,j )
2
∆t
ri ∆r2 i+ 2
and in spherical coordinates
wi,j+1 − wi,j
1
2
2
2
= 2 2 (ri+
1 wi+1,j − 2ri wi,j + r
i− 12 wi−1,j )
2
∆t
ri ∆r
where ri ; i = 1, 2, ..., N are the nodal positions along the radial direction, ri± 1 are positions
2
of points located half-way between neighboring nodes, ∆r is the mesh spacing and j is the
time step index. Note r1 = 0 is the node located at the center of the cylinder (or sphere).
All ideas presented before can be directly applied here with little modification. However,
special care is required to handle the symmetry condition at the origin r = 0. For symmetry,
it is required that
∂u
=0
∂r
A second order order accurate finite difference formula at the origin can be obtained
using a fictitious node located at r = −∆r giving
w2,j − w−2,j
=0
2∆r
i.e. w2,j = w−2,j , the condition of radial symmetry around the origin. The fictitious node
value w−2,j is then eliminated by combining the above with the finite difference formula
being used for the interior nodes giving the result
w1,j+1 − w1,j
2
=
(w2,j − w1,j )
∆t
∆r2
Alternatively and equivalently, by means of a Maclaurin expansion one can show that at
the origin, the following form of the diffusion equation is valid (when one has symmetry at
the origin),
ut = 2urr
for two-dimensional problems in cylindrical coordinates and
ut = 3urr
for three-dimensional problems in spherical coordinates. These expressions can then be
readily discretized to obtain finite difference formulae for the point at the origin.
15
If no symmetry can be assumed the following expressions can be used to approximate
the Laplacian,
∇2 u ≈
4(wM,j − w0,j )
∆r2
∇2 u ≈
6(wM,j − w0,j )
∆r2
for cylindrical systems and
for spherical systems. Here wM,j is the nearest-neighbor mean value of w obtained by averaging over all nearest neighbor nodes to the node at the origin. The above approximations
can then be used in the original PDE to obtain finite difference formulae for w0,j .
3.6
Multidimensional Problems
Multidimensional problems for parabolic PDEs result if more than one spatial variable is
involved. For instance, for the two-dimensional plate a ≤ x ≤ b, c ≤ y ≤ d, the linear
diffusion equation is
∂u
∂ 2u ∂ 2u
= α[ 2 + 2 ]
∂t
∂x
∂y
this must be solved for u(x, y, t) subject to suitable boundary and initial conditions. A
particularly simple example would require
u=0
for x = a and x = b for all y and for y = a, y = b, for all x. This starting from
u(x, y, 0) = f(x, y)
This problem can also be solved easily by separation of variables.
The finite difference method proceeds as before by first creating a mesh in space subdividing the system into M subdomains in x and N in y yielding area elements of size ∆x × ∆y
where ∆x = (b − a)/M and ∆y = (d − c)/N. Next, a mesh of spacing ∆t is also created
along the time axis. The method then computes approximations to the values of u at the
nodal locations, i.e. wi,j,k ≈ u(xi, yj , tk ).
The simples scheme is the explicit one given by
wi,j,k+1 − wi,j,k
wi+1,j,k − 2wi,j,k + wi−1,j,k wi,j+1,k − 2wi,j,k + wi,j−1,k
− α[
]=0
∆t
∆x2
∆y 2
Note that each equation in the above set uses the values of w at five contiguous nodes
i, j; i − 1, j; i + 1, j; i, j + 1; i, j − 1) for time level k to compute the value of wi,j,k+1 . As in
16
the one-dimensional case, the explicit scheme is conditionally stable and the CFL restriction
in this case requires that
∆t ≤
1 ∆x2∆y 2
2α ∆x2 + ∆y 2
Several implicit schemes are possible. For instance, if the two spatial dimensions are
handled implicitly one obtains the scheme
wi,j,k+1 − wi,j,k
wi+1,j,k+1 − 2wi,j,k+1 + wi−1,j,k+1
− α[
+
∆t
∆x2
wi,j+1,k+1 − 2wi,j,k+1 + wi,j−1,k+1
+
]=0
∆y 2
Note that each equation in the above set uses the values of w at the five contiguous nodes
i, j; i−1, j; i+1, j; i, j+1; i, j−1) for time level k+1 to compute the value of wi,j,k+1 . Although
the matrix associated with the resulting system is still sparse, at its best (depending on clever
ordering of the nodes) it is pentadiagonal.
Alternatively, one may decide to handle only one dimension implicitly (say x) during
the first part of the time step and switch the choice for the second part. Since only three
contiguous nodes at the same time level are involved at each implicit step, a much easier to
solve tridiagonal system is obtained. This is called th Alternating Direction Implicit (ADI)
scheme. Therefore, for the first half of the time step ∆t one uses
wi+1,j,k+1/2 − 2wi,j,k+1/2 + wi−1,j,k+1/2
wi,j,k+1/2 − wi,j,k
− α[
+
1
∆t
∆x2
2
wi,j+1,k − 2wi,j,k + wi,j−1,k
+
]=0
∆y 2
to calculate wi,j,k+1/2 and this is followed by
wi,j,k+1 − wi,j,k+1/2
wi+1,j,k+1/2 − 2wi,j,k+1/2 + wi−1,j,k+1/2
− α[
+
1
∆t
∆x2
2
wi,j+1,k+1 − 2wi,j,k+1 + wi,j−1,k+1
+
]=0
∆y 2
to calculate wi,j,k+1 .
Exercise. Consider the two dimensional formed by the unit square. Let f(x, y) = 1
and assume u = 0 at all four boundaries of the square. Solve the problem using the explicit
method and the alternating direction implicit method and estimate the accuracy of the
calculation in each case.
4
Hyperbolic Problems: Finite Difference Methods
Hyperbolic equations occur in many engineering and scientific applications. They are frequently associated with conservation principles, propagation problems and vibrations.
17
4.1
First Order Hyperbolic Partial Differential Equations
The most common example of first order hyperbolic PDE is the so-called conservation law
∂u ∂f(u)
+
=0
∂t
∂x
which must be solved for a given f(u), subject to the initial condition
u(x, 0) = u0 (x)
and a sometimes a suitable boundary condition in order to obtain u(x, t) . An example of
the above is what happens when motorists driving at steady speed on a highway approach
an area with a traffic jam. As autos approach the jammed area they are forced to stop and
join the jam. The size of the jam increases as the jammed area propagates in the opposite
direction of the traffic.
By making
a = a(u) =
∂f
∂u
one obtains the so-called advection equation
∂u
∂u
+a
=0
∂t
∂x
A key feature of the advection equation is that u has constant values along certain sets
of points on the x − t plane called characteristics. I.e. along such characteristic directions
du
∂u
dx ∂u
=
+( )
dt
∂t
dt ∂x
The characteristics are solutions of the following ODE
dx
= a(x, t)
dt
Therefore, given a collection of points xj at time t = 0, characteristic directions can be
determined by numerical solution of the the resulting initial value problem for x(t). For
instance, if a is simply a constant the solution of the advection equation becomes
u(x, t) = u0(x − at)
a forward moving wave.
18
4.2
Finite Difference Methods for Advection Problems
For numerical solution of the advection equation using finite difference methods one starts by
subdividing the space domain into N subdomains represented by nodes xi with i = 1, 2, ...N
and also creating a mesh in time with nodes tj , with j = 1, 2, ...M
A numerical approximation wi,j ≈ u(xi , tj ) to the solution of the advection equation can
be obtained using an explicit finite difference scheme as follows (assuming a = constant)
wi,j − wi−1,j
wi,j+1 − wi,j
+a
=0
∆t
∆x
Note that forward differences have been used for both time and space. Introducing ν =
a∆t/∆x, the scheme becomes
wi,j+1 = (1 − ν)wi,j + νwi−1,j
This shows that the computed value of wi,j+1 at time level j + 1 depends on the values of wi,j
and wi−1,j . In turn, the value of wi,j depends on the values of wi,j−1 and wi−1,j−1 and that of
wi−1,j on those of wi−1,j−1 and wi−1,j−1 and so on forming a triangular domain of dependence.
For convergence of the numerical scheme it is necessary that its domain of dependence
contains the domain of dependence of the original PDE. This is the Courant-Friedrichs-Lewy
condition for the convergence of the explicit scheme for the advection equation. Specifically,
for a > 0, the CFL condition requires that
ν=a
∆t
≤1
∆x
If seeking greater accuracy, one uses instead of a forward difference formula for the
space derivative a central difference formula, unexpectedly, an absolutely unstable scheme is
obtained!
A three point scheme with improved convergence characteristics is the upwind scheme in
which the approximation is determined by the sign of a, i.e.
wi,j+1 =
(
(1 + ν)wi,j − νwi+1,j if a < 0
(1 + ν)wi,j + νwi−1,j if a > 0
Exercise. Use the upwind scheme to solve the hyperbolic equation
ut + a(x, t)ux = 0
with
a(x, t) =
1 + x2
1 + 2xt + 2x2 + x4
and subject to
u(x, 0) = 1
19
for 0.2 ≤ x ≤ 0.4 and 0 otherwise, and
u(0, t) = 0
A more accurate method is the Lax-Wendroff scheme given by
1
1
wi,j+1 = ν(1 + ν)wi−1,j + (1 − ν 2 )wi,j − ν(1 − ν)wi+1,j
2
2
As can be shown by Fourier mode analysis, this scheme is stable for |ν| ≤ 1.
Exercise. Use the Law-Wendroff scheme to solve the problem in the previous exerciase
and compare your results.
The box method is an implicit scheme used to solve the advection equation. The scheme
is given as follows
j+1/2
j+1/2
wi+1,j+1 = wi,j + (1 + νi+1/2 )−1 (1 − νi+1/2 )(wi+1,j − wi,j+1 )
When applied in practice the scheme is used under conditions that guarantee unconditional
stability (i.e. ν > 0).
Another useful method is the leapfrog scheme given by
wi,j+1 = wi,j−1 − (a
∆t
)(wi+1,j − wi−1,j )
∆x
Note that this is an explicit scheme requiring special handling of the starting conditions (e.g.
Lax-Wendroff). Stability of this scheme also requires |ν| ≤ 1.
4.3
Finite Difference Methods for Second Order Hyperbolic Partial Differential Equations
The wave equation is the most common example of second order hyperbolic partial differential equations. The goal is to determine the function u(x, t) satisfying
2
∂ 2u
2∂ u
(x,
t)
=
c
(x, t) =
∂t2
∂x2
on 0 < x < l, 0 < t and subject to the conditions
u(0, t) = u(l, t) = 0,
t>0
u(x, 0) = f(x), 0 ≤ x ≤ l
∂u
(x, 0) = g(x) 0 ≤ x ≤ l
∂t
The implementation of the finite difference method in this case proceeds as follows. First
a mesh is set up by subdiving time and space in subintervals bounded by mesh points
20
xi = (i − 1)h; i = 1, 2, ..., N and tj = (j − 1)k; j = 1, 2, ..., M. Second order accurate finite
difference formulae for all the interior points give
wi+1,j − 2wi,j + wi−1,j
wi,j+1 − 2wi,j + wi,j−1
− c2
=0
2
k
h2
which can be solved for wi,j+1 to give
wi,j+1 = 2[1 − (
αk 2
ck
) ]wi,j + ( )2(wi+1,j + wi−1,j ) − wi,j−1
h
h
for i = 1, .., N and j = 1, 2, ....
A second order Taylor polynomial can now be constructed to satisfy the derivative initial
condition. The result is
wi,2 = wi,1 + kg(xi ) +
c2 k 2 f(xi+1 − 2f(xi ) + f(xi−1 )
[
]
2
h2
As before, the explicit method is conditionally stable. The CFL restriction in this case is
ck/h ≤ 1. Unconditionally stable implicit methods for the wave equation are also available.
5
Exercises
1.- Find approximate solutions u(x, y) to the following problem in x ∈ [0, 1], y ∈ [0, 1],
∂ 2u ∂ 2u
+
=0
∂x2 ∂y 2
subject to
u(0, y) = u(1, y) = u(x, 0) = 0
u(x, 1) = sin(πx)
using the separation of variables method.
2.- Find approximate solutions u(x, t) to the following problem in x ∈ [0, 1], t > 0,
∂ 2u ∂u
−
=0
∂x2
∂t
subject to
u(0, t) = u(1, t) = 0
u(x, 0) = 1
21
using the separation of variables method.
3.- Find approximate solutions u(x, t) to the following problem in x ∈ [0, 1], t > 0,
∂ 2u ∂u
−
=0
∂x2
∂t
subject to
u(0, t) = u(x, 0) = 0
u(0, t) = t
22
Fly UP