chap5.html

5.1 Diffusion Equations

The diffusion equation,

∂ψ

∂t

= ∇. (D ∇ψ) + s,

(5.1)

arises in heat conduction, neutron transport, particle diffusion , and numerous other situations. There is a clear difference between the time variable, t, and the spatial variables x. We'll talk mostly for brevity as if there is only one spatial dimension x, but this discussion can readily be generalized to multiple spatial dimensions. The highest time derivative is ∂/∂t, first order. The highest spatial derivative is second order ∂²/∂x². The equation is classified as parabolic.

Figure 5.1: Solving a diffusion problem generally requires a combination of spatial boundary conditions and temporal initial condition. Then the solution is propagated upward (forward in time), to fill in the multi-dimensional (time and space) solution domain.

Consequently, boundary conditions are not applied all around a closed contour (in the x-t plane) but generally only at the two ends of the spatial range, and at one "initial" time, as illustrated in Fig. 5.1. The dependent variable solution is propagated from the initial condition forward in time (conventionally drawn as the vertical direction).

5.2 Time Advance Choices and Stability

5.2.1 Forward Time, Centered Space

For simplicity we'll take D to be uniform in one cartesian dimension, and the meshes to be uniform with spacing ∆x and ∆t. Then one way to write the equation in discrete finite differences is

ψ⁽ⁿ⁺¹⁾_j −ψ⁽ⁿ⁾_j

∆t

= D

ψ⁽ⁿ⁾_j+1−2ψ⁽ⁿ⁾_j+ψ⁽ⁿ⁾_j−1

∆x²

+ s⁽ⁿ⁾_j.

(5.2)

We use (n) to denote time index, and put it as a superscript to distinguish it from the space indices j [k, l]. Of course this notation is not raising to a power. Notice in this equation the second order derivative in space is naturally centered and symmetric. However, the time derivative is not centered in time. It is really the value at n+1/2, not at the time index of everything else: n. This scheme is therefore "Forward" in Time, but Centered in Space (FTCS); see Fig. 5.2.

Figure 5.2: Forward Time, Centered Space, (FTCS) difference scheme.

We immediately know from our previous experience that, because it is not centered in time, this scheme's accuracy is going to be only first order in ∆t. Also, this scheme is explicit in time. The ψ at n+1 is obtained using only prior (n) values of the other quantities:

ψ⁽ⁿ⁺¹⁾_j = ψ⁽ⁿ⁾_j +

D∆t

∆x²

(ψ⁽ⁿ⁾_j+1−2ψ⁽ⁿ⁾_j+ψ⁽ⁿ⁾_j−1) + ∆t s⁽ⁿ⁾_j.

(5.3)

A question then arises as to whether this scheme is stable. For an ordinary differential equation, we saw that with explicit integration, there was a maximum step size that could be allowed before the scheme became unstable. The same is true for hyperbolic and parabolic partial differential equations. For stability analysis, we ignore the source S (because we are really analysing the deviation of the solution²⁷). However, even so, it's a bit difficult to see immediately how to evaluate the amplification factor, because for partial differential equations there is variation in the spatial dimension(s) that has to be accounted for. It wasn't present for ordinary differential equations. The way this is generally handled is to turn the partial differential equation into an ordinary differential equation by examining separately all the Fourier components of the spatial variation. This sort of analysis is called Von Neumann stability analysis. It gives a precisely correct answer only for uniform grids and coefficients, but it is usually approximately correct, and hence in practice very useful even for non-uniform cases.

A Fourier component varies in space like exp(i k_xx) where k_x is the wave number in the x-direction, (and i is here the square root of minus 1). For such a Fourier component, ψ_j ∝ exp(ik_x∆x j), so that ψ_j+1 = exp(ik_x∆x)ψ_j and ψ_j−1 = exp(−ik_x∆x)ψ_j. Therefore

ψ_j+1−2ψ_j+ψ_j−1 = ( e^ik_x∆x−2+e^−ik_x∆x ) ψ_j = −4sin²(

k_x∆x

)ψ_j .

(5.4)

Then substituting into eq. (5.3) for this Fourier component we find

ψ⁽ⁿ⁺¹⁾ =

⎡
⎣

1 −

D∆t

∆x²

4sin²(

k_x∆ x

)

⎤
⎦

Amplification factor

ψ⁽ⁿ⁾.

(5.5)

The amplification factor from each step to the next is the square bracket term. If it has a magnitude greater than 1, then instability will occur. If D is negative it will in fact be greater than 1. This instability is not a numerical instability, though. It is a physical instability. The diffusion coefficient must be positive otherwise the diffusion equation is unstable regardless of numerical methods. So D must be positive; and so are ∆t, ∆x. Therefore numerical instability will arise if the magnitude of the second (negative) term in the amplification factor exceeds 2.

If k_x∆x is small, then that will make the second term small and unproblematic. We are most concerned about larger k_x values that can make sin²(k_x∆x/2) approximately unity. In fact, the largest k_x value that can be represented on a finite grid²⁸ is such that the phase difference (k_x∆x) between adjacent values is π radians. That corresponds to a solution that oscillates in sign between adjacent nodes. For that Fourier component, therefore, sin²(k_x∆x/2) = 1.

Stability requires all Fourier modes to be stable, including the worst mode that has sin²(k_x∆x/2) = 1. Therefore the condition for stability is

4D∆t

∆x²

< 2.

(5.6)

There is, for the FTCS scheme, a maximum stable timestep equal to ∆x²/2D.

Incidentally, the fact that ∆t must therefore be no bigger than something proportional to ∆x² makes the first order accuracy in time less of a problem. In fact, for a timestep at the stability limit, as we decrease ∆x, improving the spatial accuracy proportional to ∆x² because of the second order accuracy in space, we also improve the temporal accuracy by the same factor, proportional to ∆x² because ∆t ∝ ∆x².

5.2.2 Backward Time, Centered Space. Implicit Scheme.

In order to counteract the instability, we learned earlier that implicit schemes are helpful. The natural implicit time advance is simply to say that we use the values of the updated variables to evaluate the rest of the equation, instead of the prior values:

ψ⁽ⁿ⁺¹⁾_j = ψ⁽ⁿ⁾_j +

D∆t

∆x²

(ψ⁽ⁿ⁺¹⁾_j+1−2ψ⁽ⁿ⁺¹⁾_j+ψ⁽ⁿ⁺¹⁾_j−1) + ∆t s⁽ⁿ⁺¹⁾_j;

(5.7)

This is a Backward in Time, Centered in Space scheme; illustrated in Fig. 5.3.

Figure 5.3: Backward time, centered space, (BTCS) difference scheme.

We'll see a little later how to actually solve this equation for the values at n+1, but we can do the same stability analysis on it without knowing. The combination for the spatial Fourier mode is just as in eq. (5.4), so the update equation (ignoring S) for a Fourier mode is

⎡
⎣

1 +

D∆t

∆x²

4sin²(

k_x∆ x

)

⎤
⎦

Inverse of Amplification factor

ψ⁽ⁿ⁺¹⁾ = ψ⁽ⁿ⁾.

(5.8)

The amplification factor is the inverse of the square bracket factor on the left. That square bracket has magnitude always greater than one. Therefore the BTCS scheme is unconditionally stable. We can take timesteps as large as we like.

It turns out, however, that the only first order accuracy of this scheme, like the FTCS scheme, means that we don't generally want much to exceed the previously calculated stability limit. If we do so, we increasingly sacrifice temporal accuracy, even though not stability.

5.2.3 Partially implicit, Crank-Nicholson schemes.

The best choice, for optimizing the efficiency of the numerical solution of diffusive problems is to use a scheme that is part forward and part backward. A combination of forward and backward, in which θ is the weight of the implicit or backward proportion, is to write

ψ⁽ⁿ⁺¹⁾_j = ψ⁽ⁿ⁾_j

⎡
⎣

D∆t

∆x²

(ψ⁽ⁿ⁺¹⁾_j+1−2ψ⁽ⁿ⁺¹⁾_j+ψ⁽ⁿ⁺¹⁾_j−1) + ∆t s⁽ⁿ⁺¹⁾_j

⎤
⎦

(1−θ)

⎡
⎣

D∆t

∆x²

(ψ⁽ⁿ⁾_j+1−2ψ⁽ⁿ⁾_j+ψ⁽ⁿ⁾_j−1) + ∆t s⁽ⁿ⁾_j

⎤
⎦

(5.9)

This is sometimes called the "θ-implicit" scheme. The amplification factor is straightforwardly

A =

1 − (1−θ)

4D∆t

∆x²

sin²(

k_x∆x

)

1 + θ

4D∆t

∆x²

sin²(

k_x∆x

)

(5.10)

If θ ≥ ¹/₂ then |A| ≤ 1, and the scheme is always stable. If θ < ¹/₂, then |A| ≤ 1 requires the stability criterion

∆t <

∆x²

2 D(1−2θ)

(5.11)

Thus the minimum degree of implicitness that guarantees stability for all sizes of timestep is θ = ¹/₂. This choice is called the "Crank-Nicholson" scheme.

It has a major advantage beyond stability. It is centered in time. That means it is second order accurate in time (as well as space). This accuracy makes it useful to take bigger steps than would be allowed by the (explicit advance) stability limit.

5.3 Implicit Advancing Matrix Method

An implicit or partially implicit scheme for advancing a parabolic equation generally results in equations containing more than one spatial point value at the updated time, for example ψ⁽ⁿ⁺¹⁾_j, ψ⁽ⁿ⁺¹⁾_j−1, ψ⁽ⁿ⁺¹⁾_j+1. Such an equation for all spatial positions can be written as a matrix equation. Gathering together the terms at n and at n+1 from eq. (5.9), it can be written

⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝

^··_·

B₋

B₀

B₊

^··_·

⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠

⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝

ψ⁽ⁿ⁺¹⁾₁

ψ⁽ⁿ⁺¹⁾_j

ψ⁽ⁿ⁺¹⁾_J

⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠

⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝

^··_·

C₋

C₀

C₊

^··_·

⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠

⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝

ψ⁽ⁿ⁾₁

ψ⁽ⁿ⁾_j

ψ⁽ⁿ⁾_J

⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠

⎛
⎜
⎜
⎜
⎜
⎜
⎜
⎝

s₁

s_j

s_J

⎞
⎟
⎟
⎟
⎟
⎟
⎟
⎠

(5.12)

or symbolically

Bψ_n+1 = Cψ_n + s,

(5.13)

where J is the total length of the spatial mesh (the maximum of j), and the coefficients are

B₀=1+ 2

D∆t

∆x²

θ, B₊=B₋=−

D∆t

∆x²

θ,

(5.14)

C₀=1−2

D∆t

∆x²

(1−θ), C₊=C₋=+

D∆t

∆x²

(1−θ),

(5.15)

and s_j = ∆t [θs⁽ⁿ⁺¹⁾_j + (1−θ)s⁽ⁿ⁾_j].

(5.16)

[Notice the ∆t scaling factor in s_j.]

We are here assuming that the source terms do not depend upon ψ. Eq. (5.13) can be solved formally by inverting the matrix B:

ψ_n+1 = B⁻¹Cψ_n + B⁻¹s.

(5.17)

Therefore the additional work involved in using an implicit scheme is that we have (effectively) to invert the matrix B, and multiply by the inverse at each timestep.

Provided the mesh is not too large, this can be a managable approach. In a situation where B and C do not change with time, the inversion²⁹ need only be done once; and each step in time involves only a matrix multiplication by B⁻¹C of the values from the previous timestep, ψ_n.

If the matrix has the tridiagonal form of eq. (5.12), where the entries are non-zero only on the diagonal and the immediately adjacent subdiagonals, then it will certainly be more computationally efficient to solve for the updated variable by elimination³⁰ rather than inverting the matrix and multiplying. This is a reflection of the fact that the matrix B is sparse. All but a very small fraction of its coefficients are zero. The fundamental problem is that the inverse of a sparse matrix is generally not sparse. Consequently, even though multiplication by the original sparse matrix actually requires only a few individual arithmetic operations, and easy short-cuts can be implemented, there are no obvious short-cuts for doing the matrix multiplication by the inverse. As long as we are implementing solutions in a mathematical system like Octave using the built-in facilities, we won't notice any difference because we are not using short-cuts.

5.4 Multiple Space Dimensions

Figure 5.4: Example solution, at a particular time, of a diffusion equation in two space dimensions. The value of ψ is visualized in this perspective plot as the height of the surface.

When there is more than one spatial coordinate dimension, as illustrated in Fig. 5.4, nothing changes formally about our method of solution of parabolic equations. What changes, however, is that we need a systematic way to turn the entire spatial mesh in multiple dimensions into a column vector like those in eq. (5.12). In other words, we must index all the spatial positions with a single index j. But generally if we have multiple dimensions, the natural (physical) indexing of the elements of the mesh is by a number of multiple indices equal to the number of dimensions: e.g. ψ_k,l,m, where k, l, m correspond to the coordinates x, y, z.

Reordering the mesh elements is not a particularly difficult task algebraically, but it requires an intuitively tricky geometrical mental shift. In general, if we have a quantity ψ_k,l,... indexed on a multidimensional mesh whose lengths in the different mesh dimensions are K, L, ... then we re-index them into a single index j as follows. Start with all indices equal to 1. Count through the first index, k = 1,...,K. Then increment the second index l and repeat k = 1,...,K. Continue this process for l=1,...,L. Then increment the next index (if any), and continue till all indices are exhausted. Incidentally, this is precisely the order in which the elements are stored in computer memory when using a language like Fortran when using the entire allocated array.

Figure 5.5: Conversion of a two-dimensional array into a single column vector using column order is like stacking all the columns below one another.

The result is that for the giant column vector ψ_j, the first K elements are the original multidimensionally indexed elements ψ_k,1,..., k=1,...,K; the next K elements, j=K+1,...,2K, are the original ψ_k,2,... and so on.³¹

In multiple dimensions, the second derivative operator (something like ∇²) is represented by a stencil, as illustrated in eq. (4.23). The importance of the reordering of the elements into a column vector is that although the components of the stencil that are adjacent in the first index (k) remain adjacent in j, those in the other index(es) (l) do not. For example the elements ψ_k,l and ψ_k,l+1 have j-indices j=K.(l−1)+k and j=K.l+k respectively; they are a distance K apart. This fact means that for multiple spatial dimensions the matrices B and C are no longer simply tridiagonal. Instead, they have an additional non-zero diagonal, a distance K from the main diagonal. If the boundary conditions are ψ = 0, each matrix is "Block-Tridiagonal", having a form like this:

L blocks

⎧
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩

⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣

⎫
⎪
⎪
⎬
⎪
⎪
⎭

Block

size

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦

⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭

J=K×L rows.

(5.18)

Here, c denotes the coefficient of the center of the stencil, and a denotes the coefficients of the adjacent points in the stencil. One can think of the total matrix as being composed of L×L blocks, each of size K×K. All the blocks are zero except the tridiagonal ones. And each non-zero block is itself tridiagonal (or diagonal). If there are further dimensions M, then the giant matrix is a tridiagonal composite of M×M blocks, each of a two dimensional (L×L block) type of eq. (5.18). And so on.

It is very important to think carefully about the boundary conditions. These occur at the boundaries of each block. Notice how the corner (boundary) entries of the extra subdiagonal blocks make zero certain coefficients of the subdiagonals of the total matrix. One important consequence of the block-tridiagonal form (5.18) is that it is not so easy to do tridiagonal elimination rather than matrix inversion.

5.5 Estimating Computational Cost

The computational cost of executing different numerical algorithms often has to be considered for large scale computing. This cost is generally expressed in the number of floating point operations (FLOPs). Roughly, this is the number of multiplications (or divisions) required. Additions are cheap.

For an N×N matrix, multiplying into a length-N column vector costs (simplemindedly) N rows times N operations per row: N² (FLOPs). By extending this to N columns, we see that multiplying two matrices costs N³. Although the process of inversion³² of a non-singular matrix, seems far more difficult to do, because the algorithm is much more complicated, it also costs roughly N³ operations . These estimates are accurate only to within a factor of two or so. But that is enough for most purposes.

Inverting or multiplying random matrices using Octave on my laptop, for N=1000 takes about 1 second. That seems amazingly fast to me, because it corresponds to about 1 FLOP per nanosecond. But that's about what can be achieved these days provided the cache contains enough of the data.

Trouble is that if we are dealing with two-dimensions in space, each of length K=100, then the side of the matrix we have to invert is J=K²=10000. A multiplication or inversion would take at least J³/10⁹ = 1000 seconds; that's a quarter of an hour. Waiting that long is frustrating, and if many inversions are needed, time rapidly runs out.

Inverting the matrix representing a two-dimensional problem takes K⁶ operations. And for a three-dimensional problem it's K⁹. This rapid escalation of the cost means that one doesn't generally approach multidimensional problems in terms of simplistic matrix inversion. In particular, it is difficult to implement an implicit scheme for advancing the diffusion equation. And it is probably more realistic just to use an explicit scheme, recognizing and observing the stability limit, eq. (5.6), on maximum time-step size.

Worked Example: Crank Nicholson Matrices

Express the following parabolic partial differential equation in two space variables (r and ϕ) and one time variable (t) as a matrix equation using the Crank-Nicholson scheme on uniform grids.

∇²ψ =

r²

∂

∂r

⎛
⎝

r²

∂ψ

∂r

⎞
⎠

r²

∂²ψ

∂ϕ²

∂ψ

∂t

(5.19)

Find the coefficients and layout of the required matrices if the boundary conditions are

∂ψ

∂r

at r=0, ψ = 0 at r=R, and periodicity in ϕ.

Let k and l be indices in the r and ϕ coordinates, and let ∆r and ∆ϕ denote their grid spacings. The finite difference form of the spatial derivative terms is

r²

∂

∂r

⎛
⎝

r²

∂ ψ

∂r

⎞
⎠

k,l

r_k²

⎛
⎝

r_k+1/2²

ψ_k+1,l−ψ_k,l

r_k+1−r_k

−r_k−1/2²

ψ_k,l−ψ_k−1,l

r_k−r_k−1

⎞
⎠

r_k+1/2−r_k−1/2

⎡
⎢
⎣

⎛
⎜
⎝

r_k+

∆r

r_k

⎞
⎟
⎠

(ψ_k+1,l−ψ_k,l)−

⎛
⎜
⎝

r_k−

∆r

r_k

⎞
⎟
⎠

(ψ_k,l−ψ_k−1,l)

⎤
⎥
⎦

∆r²

⎡
⎣

⎛
⎝

∆r

2r_k

⎞
⎠

ψ_k+1,l−

⎛
⎝

∆r²

2r_k²

⎞
⎠

ψ_k,l+

⎛
⎝

1−

∆r

2r_k

⎞
⎠

ψ_k−1,l

⎤
⎦

∆r²

r²

∂²ψ

∂θ²

⎢
⎢

k,l

r_k²

⎛
⎝

ψ_k,l+1−ψ_k,l

θ_l+1−θ_l

−

ψ_k,l−ψ_k,l−1

θ_l−θ_l−1

⎞
⎠

θ_l+1/2−θ_l−1/2

(ψ_k,l+1−2ψ_k,l+ψ_k,l−1)

r_k²∆θ²

(5.20)

Therefore, expressing the positions across the spatial grid in terms of a single index j=k+K.(l−1) (where K is the size of the k-grid, and L of the l-grid), the differential operator ∇² becomes a matrix M multiplying the vector ψ =ψ_j, of the form

M =

⎡
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣

⎫
⎪
⎪
⎬
⎪
⎪
⎭

Block

size

⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦

⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎭

L blocks.

(5.21)

The explanation of this form is as follows. The coefficients of a generic row corresponding to r-position index k are

⎛
⎝

1−

∆r

2r_k

⎞
⎠

∆r²

, b=

⎛
⎝

∆r

2r_k

⎞
⎠

∆r²

, d=

r_k²∆ϕ²

, c=−(a+b+2d).

(5.22)

The periodic boundary conditions in ϕ are implemented by the appearance of off-diagonal, d-type, blocks at the upper right and lower left of the matrix. The r-boundary conditions are the same for all the blocks on the diagonal, the tridiagonal acb-type blocks. At the r=R boundary, (which would be the bottom row of each block), the condition ψ = 0 means no contribution arises to the differential operator from the ψ-value there. The condition therefore allows us simply to omit the r=R row from the matrix, choosing the k=K index to refer to the position r=R−∆r. At the r=0 end (the top row of each block), the condition

∂ψ

∂ r

can be implemented in a properly centered way by choosing the r-grid to be aligned to the half-integral positions.³³ In other words, r₀=−∆r/2, r₁=+∆r/2, r_n=(n−1/2)∆r. (The r₀ position values are not represented in the matrices.) In that case, there is zero contribution to the difference scheme from the first derivative at position r_1−1/2=0 (because r² there is zero), and the r-second-derivative operator at k=1 becomes

(1+

∆r

2r₁

)²(ψ_2,l−ψ_1,l)/∆r²=4(ψ_2,l−ψ_1,l)/∆r²

, giving rise to equal and opposite coefficients e=4/∆r². The diagonal entry on those k=1 rows is minus the sum of all the other coefficients on the row: f=−(e+2d).

The Crank-Nicholson scheme for the differential equation time advance is then

ψ⁽ⁿ⁺¹⁾−ψ⁽ⁿ⁾=

∆t

Mψ⁽ⁿ⁺¹⁾ +

∆t

Mψ⁽ⁿ⁾,

(5.23)

which on rearrangement becomes

ψ⁽ⁿ⁺¹⁾ = (I−

∆t

M)⁻¹ (I+

∆t

M)ψ⁽ⁿ⁾.

(5.24)

Exercise 5. Diffusion and Parabolic Equations.

1. Write a computer code to solve the diffusive equation

∂ψ

∂t

= D

∂² ψ

∂ x²

+ s(x)

For constant, uniform diffusivity D and constant specified source s(x). Use a uniform x-mesh with N_x nodes. Consider boundary conditions to be ψ fixed, equal to ψ₁, ψ₂ at the domain boundaries, x=−1,1, and the initial condition to be ψ = 0 at t=0.

Construct a matrix G=G_ij such that Gψ = ∇²ψ. Use it to implement the FTCS scheme

ψ⁽ⁿ⁺¹⁾ = (I+∆t G)ψ⁽ⁿ⁾ + ∆t s,

paying special attention to the boundary conditions.

Solve the time-dependent problem for D=1, s=1, N_x=30, ψ₁=ψ₂=0, t=0→1, storing your results in a matrix ψ(x,t)=ψ_{j_x,j_t}, and plotting that matrix at the end of the solution, for examination.

Experiment with various ∆t to establish the dependence of the accuracy and stability of your solution on ∆t. In particular, without finding an "exact" solution of the equation by any other method,

(i) Find experimentally the value of ∆t above which the scheme becomes unstable.

(ii) Estimate experimentally the fractional error arising from finite time step duration in ψ(t=1,x=0) when using a ∆t approximately equal to the maximum stable value.

(iii) By varying N_x, estimate experimentally the fractional error at N_x=30 arising from finite spatial differences. Which is more significant, time or space difference error?

2. Develop a modified version of your code to implement the θ-implicit scheme:

(I−∆tθG) ψ⁽ⁿ⁺¹⁾ = (I+∆t (1−θ) G)ψ⁽ⁿ⁾ + ∆t s,

in the form

ψ⁽ⁿ⁺¹⁾ = B⁻¹C ψ⁽ⁿ⁾ + B⁻¹ ∆t s

(i) Experiment with different ∆t and θ values, for the same time-dependent problem and find experimentally the value of θ for which instability disappears for all ∆t.

(ii) Suppose we are limited to only 50 timesteps to solve over the time 0 < t ≤ 1 so ∆t=0.02. Find experimentally the optimum value of θ which produces the most accurate results.

Submit the following as your solution for each part:

Your code in a computer format that is capable of being executed, citing the language it is written in.
The requested experimental ∆t and/or θ values.
A plot of your solution for at least one of the cases.
A brief description of how you determined the accuracy of the result.

HEAD

Chapter 5 Diffusion. Parabolic Partial Differential Equations

5.1 Diffusion Equations

5.2 Time Advance Choices and Stability

5.2.1 Forward Time, Centered Space

5.2.2 Backward Time, Centered Space. Implicit Scheme.

5.2.3 Partially implicit, Crank-Nicholson schemes.

5.3 Implicit Advancing Matrix Method

5.4 Multiple Space Dimensions

5.5 Estimating Computational Cost

Worked Example: Crank Nicholson Matrices

Exercise 5. Diffusion and Parabolic Equations.

Chapter 5
Diffusion. Parabolic Partial Differential Equations