6 Numerical Solution of Parabolic Equations

(1)

6 Numerical Solution of Parabolic Equations

6.1 Some Theory and Examples

Our first example is the so-called initial value problem for the heat equation (u = temperature in an infinitely long wire):

u

t

= u

xx

in Ω = { (x, t) : x ∈

R

, 0 < t < T } , u(x, 0) = f(x).

We define the

singularity function

S(x, t) := 1

√ 4πt exp − x

²

4t

which solves at least the differential equation for t 6= 0.

Parabolic Problems TU Bergakademie Freiberg 211

Plot of the singularity function

−1

−0.5 0

0.5 1

0 0.02 0.04 0.06 0.08 0.1

0 1 2 3 4

t x

S(x,t)

−1 −0.5 0 0.5 1

0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

x

t 0.5

0.5 0.5

0.5

1

11

1

1.5

3

Remark: For n spatial coordinates the singularity function is S(x, t) = 1

(4πt)

^n/2

exp( −kxk

²2

4t ).

If f satisfies the growth condition

| f(x) | ≤ M exp(Ax

²

) with some constants M ≥ 0, A > 0, so

u(x, t) =

(R_∞

−∞

S(x − y, t)f(y) dy, 0 < t < T,

f (x), t = 0,

solves the IVP, if T <

_A⁴

(Poisson’s formula). In particular, u ∈ C( ¯ Ω) and u

t

, u

xx

∈ C(Ω).

The IVP is well-posed.

(2)

As a second example we consider the parabolic initial-boundary value problem,

u

t

= u

xx

for t > 0, 0 < x < 1, u(0, t) = u(1, t) = 0 for t > 0,

u(x, 0) = u

0

(x) for 0 ≤ x ≤ 1,

which modelles the heat flow in a homogeneous thin rod with fixed zero temperature at both ends.

If the initial condition has can be described as a Fourier series in the form

u

0

(x) =

X∞ j=1

α

j

sin(jπx), (1)

then the solution of the IBVP is u(x, t) =

X∞ j=1

α

j

exp( − (jπ)

²

t) sin(jπx).

Since the functions {sin(jπx)}

^∞j=1

form a complete orthogonal system of the function space L

²

(0, 1), each function u

0

∈ L

²

(0, 1) can be written in the form (1).

The coefficients are given by α

j

= 2

Z1 0

u

0

(x) sin(jπx) dx.

Example: For u

0

(x) = 1 − 2

x

−

¹₂

the coefficients are given by α

j

=

_j2⁸π²

sin

^jπ₂

, (j ∈

N

). Visualisation of some partial sums:

0 0.2 0.4 0.6 0.8 1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

N=1 N=3 N=5 N=19 u0(x)

0 0.5 1

0 0.2 0.4 0.6 0.8 1 0

0.2 0.4 0.6 0.8 1

x t ⁰ ¹ ² ³ ⁴ ⁵ ⁶ ⁷

10⁻¹²⁰ 10⁻¹⁰⁰ 10⁻⁸⁰ 10⁻⁶⁰ 10⁻⁴⁰ 10⁻²⁰ 10⁰

m

t=0 t=1/2

Fourier series solution of IBVP (2) u

0

(x) = 1 − 2 | x −

¹₂

| in domain (x, t) ∈ [0, 1] × [0, 1] (left).

Decay of first 4 nonzero Fourier coefficients of u(t) (right).

(3)

Maximum principle for the heat equation

The uniqueness of the solution and the continuous dependence of the solution the data are a consequence of the following maximum principle:

Let Ω := (0, 1) × (0, T ) and u ∈ C( ¯ Ω) with u

t

, u

xx

∈ C(Ω) be a solution of the heat equation u

t

= u

xx

in Ω.

Then, the maximum of u can be found on the boundary of Ω, i. e. for x = 0, x = 1 or t = 0.

Wy does this maximum principle indeed ensure, that our IBVP is well posed?

Neumann boundary conditions

If we replace the homogeneous Dirichlet boundary conditions by homogeneous Neumann conditions

u

x

(0) = u

x

(1) = 0, and if the initial condition can be written in the form

u

0

(x) =

X∞ j=0

α

j

cos(jπx), the solution is given by

u(x, t) =

X∞ j=0

α

j

exp(−(jπ)

²

t) cos(jπx).

6.2 Numerics of the one-dimensional model problem In the following we consider the initial boundary value problem (IBVP) modelling heat flow in a thin rod with a given time dependent temperature at the endpoints:

u

t

= κu

xx

, x ∈ (0, 1), t > 0, (2a)

u(0, t) = g

0

(t), t > 0, (2b)

u(1, t) = g

1

(t), t > 0, (2c)

u(x, 0) = u

0

(x), x ∈ [0, 1]. (2d)

The constant κ is the heat conductivity (which we set to one in the following).

Remark: In 2D or 3D the heat equation would be u

t

= κ∆

x

u. One

would set initial data (for t = 0) and boundary data (on ∂Ω).

(4)

We first discretize the IBVP (2) in the spatial variable x only, leaving the time t continuous.

To this end we proceed as in the elliptic case and introduce the grid points

0 = x

0

< x

1

< · · · < x

J

< x

J+1

= 1 using a fixed grid spacing ∆x = 1/(J + 1) and approximate

u

xx

|

x=xj

≈ [A

∆xu]j

, j = 1, 2, . . . , J, with

A

∆x

= 1

∆x

²

tridiag(1, − 2, 1). (3)

If

u

=

u(t)

denotes the vector with components

u

j

(t) ≈ u(x

j

, t), t > 0, j = 1, 2, . . . , J, then (2) is transformed into the semi-discrete system of ODEs

u⁰

(t) = A

∆xu(t) +g(t),

(4a)

u(0) =u0

, (4b)

with [u

0

]

j

= u

0

(x

j

), j = 1, 2, . . . , J as well as

g(t) =

1 ∆x

²

[g

0

(t), 0, . . . , 0, g

1

(t)]

^T

∈

R^J

.

Obviously, (4) is an IVP for an inhomogeneous linear system of ODEs with constant coefficient matrix A

∆x

.

We can now solve (4) with known numerical methods for solving ODEs.

Introducing the fixed time step ∆t > 0, we set

U

_jⁿ

≈ [u(t

n

)]

j

≈ u(x

j

, t

n

), t

n

= n∆t.

The approximation of the solution of a time-dependent PDE as a system of ODEs along the “lines” { (x

j

, t) : t > 0 } is known as the method of lines.

We will consider in particular the explicit and implicit Euler method and

the trapezoidal rule for solving the resulting ODE system.

(5)

Applying the explicit Euler method to (4) (setting g

0

(t) = g

1

(t) = 0 for now) leads to

U

_jⁿ⁺¹

= U

_jⁿ

+ ∆t

∆x

²

U

_j−1ⁿ

− 2U

_jⁿ

+ U

_j+1ⁿ

, 1 ≤ j ≤ J, n = 0, 1, 2 . . . . This corresponds to the finite difference approximation

U

_jⁿ⁺¹

− U

_jⁿ

|

∆t

{z }

≈ut

= U

_jⁿ₋₁

− 2U

_jⁿ

+ U

_j+1ⁿ

∆x

²

| {z }

≈uxx

(5)

of the differential equation (2a).

In matrix notation:

Uⁿ⁺¹

= (I + ∆t A

∆x

)U

ⁿ

, n = 0, 1, 2, . . . , with

Uⁿ

= [U

₁ⁿ

, U

₂ⁿ

, . . . , U

_Jⁿ

]

^>

.

As usual, we define the local discretization error of the difference scheme (5) to be the residual obtained on inserting the exact solution into the difference scheme:

d(x, t) := u(x, t+∆t) − u(x, t)

∆t − u(x − ∆x, t) − 2u(x, t) + u(x+∆x, t)

∆x

²

.

Using Taylor expansions in (x, t) one easily obtains:

d

^eE

(x, t) = ∆t

2 u

tt

− ∆x

²

12 u

xxxx

+ O(∆t

²

) + O(∆x

⁴

). (6) Terminology: the explicit Euler method for solving the heat equation is consistent of first order in time and of second order in space.

Besides the asymptotic statement (6) we also require upper bounds for d

^eE

(x, t). Truncating the Taylor expansions with a remainder term, we obtain in time

u(x, t + ∆t) = (u + ∆t u

t

) |

(x,t)

+ ∆t

²

2 u

tt

(x, τ), τ ∈ (t, t + ∆t).

Proceeding analogously in x yields d

^eE

(x, t) = ∆t

2 u

tt

(x, τ) − ∆x

²

12 u

xxxx

(ξ, t), ξ ∈ (x − ∆x, x + ∆x), and we obtain, setting µ :=

_∆x^∆t2

,

|d

^eE

(x, t)| ≤ ∆t

2 M

tt

+ ∆x

²

12 M

xxxx

= ∆t 2

M

tt

+ 1 6µ M

xxxx

, (7)

assuming |u | ≤ M und |u | ≤ M on [0, 1] × [0, T ].

(6)

Applying the implicit Euler method to (4) one obtains, in place of (5), the implicit difference scheme

U

_jⁿ⁺¹

− U

_jⁿ

∆t = U

_jⁿ⁺¹₋₁

− 2U

_jⁿ⁺¹

+ U

_j+1ⁿ⁺¹

∆x

²

. (8)

The calculation of

Uⁿ⁺¹

from

Uⁿ

is seen to require the solution of a linear system of equations with coefficient matrix I − ∆tA

∆x

. Here Taylor expansion results in a local discretization error of

d

^iE

(x, t) = − ∆t

2 u

t

+ ∆x

²

12 u

xxxx

+ O(∆t

²

) + O(∆x

⁴

).

Applying the trapezoidal rule, which for an ODE

y⁰

(t) =

f(t,y(t))

is given by

yⁿ⁺¹

=

yⁿ

+ ∆t 2

f(tn

,

yⁿ

) +

f(tn+1

,

yⁿ⁺¹

) yields another implicit scheme

U

_jⁿ⁺¹

− U

_jⁿ

∆t = 1 2

U

_jⁿ₋₁

− 2U

_jⁿ

+ U

_j+1ⁿ

∆x

²

+ U

_jⁿ⁺¹₋₁

− 2U

_jⁿ⁺¹

+ U

_j+1ⁿ⁺¹

∆x

²

, which in this context is known as the Crank-Nicolson

¹

scheme.

This method requires in each time step the solution of a linear system of equations with the coefficient matrix I −

^∆t₂

A

∆x

.

For Crank-Nicolson (CN) there holds

d

^CN

(x, t) = O(∆t

²

) + O(∆x

²

).

1J. Crank and P. Nicolson, 1947

6.3 Stability and Convergence 6.3.1 Absolute Stability

In the methods of lines, an ODE solver is applied to the semidiscrete system (4).

Recalling the regions of absolute stability of the explicit and implicit Euler and trapezoidal methods, we are led to investigate the behavior of the eigenvalues of the matrix A

∆x

in (3) as ∆x → 0.

−3 −2 −1 0 1

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2

explicit Euler R(ˆλ) = 1 + ˆλ

−1 0 1 2 3

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2

implicit Euler R(ˆλ) = ¹

1−λˆ

−2 −1 0 1 2

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2

trapezoidal R(ˆλ) =^1+ˆ^λ/2

1−λ/2ˆ

(7)

The eigenvalues of A

∆x

(cf. pg. 200) are given by λ

j

= − 4

∆x

²

sin

²

jπ∆x

2 , j = 1, 2 . . . , J.

The eigenvalue of largest magnitude is λ

J

= − 4

∆x

²

sin

²

π

2 (1 − ∆x)

= − 4

∆x

²

+ π

²

+ O ∆x

²

. Therefore, in order to satisfy

| R ∆t λ

j

| ≤ 1 for all j = 1, 2, . . . , J in the explicit Euler method, it is necessary that

∆t

∆x

²

≤ 1

2 . (9)

Both implicit Euler and the trapezoidal rule are absolutely stable.

Therefore, as all eigenvalues of A

∆x

lie in the left half plane, the requirement of absolute stability places no constraints on ∆t.

6.3.2 The Lax-Richtmyer equivalence theorem

All three methods considered so far are consistent, i. e., at all points (x, t) of the domain we have d(x, t) → 0 as ∆x → 0 and ∆t → 0.

To analyze their convergence, we proceed as in the case of numerical methods for ODEs and consider a finite time interval t ∈ [0, T ], T > 0 as well as a sequence of grids with grid spacings ∆x → 0, ∆t → 0.

We determine whether at every fixed grid point (x

j

, t

n

) also the global error u(x

j

, t

n

) − U

_jⁿ

tends to zero uniformly.

Problem: A sequence of grid spacings { (∆x)

k

, (∆t)

k

} can approach the point (0, 0) in different ways. Will any „refinement curve“ in the (∆x, ∆t)-plane lead to a convergent method?

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

∆ x

∆ t

∆ t = ∆ x²/4

∆ t = ∆ x²/2

∆ t = ∆ x²

∆ t = ∆ x

∆ t = ∆ x/2

Different refinement curves ∆t,∆x → 0:

red: ∆t/∆x

²

= const, blue: ∆t/∆x = const.

(8)

The appropriate stability concept for the convergence analysis of numerical methods for solving IVPs for PDEs was developed by Lax

²

and Richtmyer

³

. We assume that ∆t is chosen as a fixed function of ∆x.

Definition 6.1.

A linear finite difference scheme of the form

Uⁿ⁺¹

= B(∆t)U

ⁿ

+

gⁿ

(∆t), n = 0, 1, . . . (10) is called Lax-Richtmyer stable, if for any fixed stopping time T there exists a constant K

T

> 0 such that

kB(∆t)

ⁿ

k ≤ K

T

, ∀∆t > 0 and n ∈

N, n∆t

≤ T.

2Peter D. Lax,∗1926

3Robert Davis Richtmyer, 1910–2003

Theorem 6.2 (Lax-Richtmyer equivalence theorem).

A consistent linear finite difference scheme (10) is convergent if, and only if, it is Lax-Richtmyer stable.

Examples:

For the explicit Euler method applied to (2) the resulting matrix B

^eE

(∆t) = I + ∆t A

∆x

is symmetric, implying

k B

^eE

(∆t) k

2

= max {| λ | : λ eigenvalue of A } . If condition (9), i. e.

∆t

∆x

²

≤ 1

is satisfied, we have kB

^eE

(∆t)k

2

≤ 1, implying that the explicit Euler 2 method is Lax-Richtmyer stable and therefore convergent.

The implicit Euler and Crank-Nicolson schemes are also Lax-Richtmyer stable but without constraints on the mesh ratio or time step. The matrices B(∆t) are given in these cases by

B

^iE

(∆t) = (I − ∆t A

∆x

)

⁻¹

, B

^CN

(∆t) =

I − ∆t

2 A

∆x

₋1

I + ∆t

2 A

∆x

.

However, since the implicit Euler scheme is consistent of order 2 in x and of order 1 in t, this suggests a time step ∆t ≈ ∆x

²

.

In contrast, the CN scheme (trapezoidal rule) is consistent of order 2 in

both x and t, which suggests a time step ∆t ≈ ∆x (much better for

small ∆x).

(9)

0 0.2 0.4 0.6 0.8 1 10⁻⁷

10⁻⁶ 10⁻⁵ 10⁻⁴ 10⁻³ 10⁻²

t

En

expl. Euler dx=1/10 expl. Euler dx=1/20 Crank Nicolson, dx = 1/10 Crank Nicolson, dx = 1/20

Error E(t

n

) = max

j

|u

ⁿj

− U

_jⁿ

| of numerical solution (2) with

g

0

= g

1

= 0 and u

0

(x) = x(1 − x) obtained with explicit Euler method with ∆t = ∆x

²

/2 and the Crank-Nicolson method with ∆t = ∆x/2.

Write a Matlab program for solving the parabolic initial-boundary value problem

u

t

= u

xx

for t > 0, 0 < x < 1, u(0, t) = u(1, t) = 0 for t > 0,

u(x, 0) = 3 sin(2πx) for 0 ≤ x ≤ 1.

with the explicit and implicit Euler and the Crank-Nicolson scheme in the time interval [0, 0.1].

Analyze the behaviour of the global discretization error for different refinement curves as ∆t = c∆x

²

(c =

¹₄

,

¹₂

, 1), ∆t = ∆x etc.

Remember that the true solution is given by (cf. pg. 214) u(x, t) = 3e

⁻^4π²^t

sin(2πx).

Remark

In all three schemes considered so far we were able to show the stronger statement k B(∆t) k ≤ 1, which is sometimes called strong stability.

However, strong stability is not necessary for Lax-Richtmyer stability.

More precisely, for Lax-Richtmyer stability it is necessary that there exists a constant α such that

kB(∆t)k ≤ 1 + α∆t (12) for all sufficiently small time steps ∆t, which is evident from the inequality

k B(∆t)

ⁿ

k ≤ (1 + α∆t)

ⁿ

≤ e

^αT

, ∀ n∆t ≤ T.

(10)

Goals of Chapter 6

You should know some typical parabolic 1D-examples and their analytic solutions.

You should be able to discretize an IBVP of type (2) in space direction and to formulate the associated ODE system.

You should be able to solve this system with the explicit/implicit Euler and the Crank-Nicolson scheme.

You should know the consistency orders in ∆t and ∆x for all three schemes and the resulting reasonable refinement curves.