Exact Solutions by Separation of Variables

For this linear problem, both the PDE and our Finite Difference approximation may be solved exactly by a method known as the separation of variables. We consider the problem

\begin{matrix}u_{t} & = u_{xx}  in 0 < x < 1, t > 0 \\ with u(0, t) & = u(1, t) = 0, u(x, 0) = u_{0}(x) \end{matrix} \}  (4.3)

which is a slight modification of (3.1). The PDE has separable solutions of the form u(x, t) = X(x)T (t) provided

XT′ = X′′T  or \frac{T′}{T} = \frac{X′′}{X} = −ω^{2},  say.

As T′/T is a function of t only, while X′′/X is a function of x only, both functions must be a constant, which we take to be negative. Then the functions X(x) and T (t) take the forms

X = A cos ωx + B sin ωx, and T = Ce^{−ω^{2}t} .

If we require X to obey the boundary conditions in (4.3), namely X(0) = X(1) = 0, we obtain non-zero solutions only if A = 0 and ω = mπ, for some integer m, so that

u = B_{m} sin mπx e^{−m^{2}π^{2}t},

for some constant B_{m}. As (4.3) is a linear problem, we may combine solutions to obtain a more general solution in the form

u(x, t) = ∑_{m=1}^{∞} B_{m} sin mπx e^{−m^{2}π^{2}t} .  (4.4)

The initial condition will be satisfied if

u(x, 0) = ∑_{m=1}^{∞} B_{m} sin mπx = u_{0}(x) .  (4.5)

Thus all we need do to obtain the solution of (4.3) is to expand the initial condition u = u_{0}(x) in a Fourier series, and substitute the appropriate values of the constants B_{m} into (4.4).

Similar methods may be used to solve the difference equation (4.1), repeated here:

U_{n}^{j+1} − U_{n}^{j} = r [θδ^{2}U_{n}^{j+1} + (1 − θ)δ^{2}U_{n}^{j}]  (4.1)

We seek solutions of the form U_{n}^{j} = X_{n}T^{j}. Substituting into (4.1) and separating the terms depending on j and n leads to

\frac{T^{j+1} − T^{j}}{r [θT^{j+1} + (1 − θ)T^{j}]} = \frac{X_{n+1} − 2X_{n} + X_{n−1}}{X_{n}} = σ,  say  (4.6)

where by similar arguments to those used above σ must be a constant. X_{n} therefore obeys the second order difference equation

X_{n+1} − (σ + 2)X_{n} + X_{n−1} = 0 ,

which (c.f. (1.8)) is solved by the substitution X_{n} = (p)^{n} if p^{2} − (σ + 2)p + 1 = 0 to give

X_{n} = A(p_{+})^{n} + B(p_{−})^{n} , where p_{±} = \frac{(σ + 2) ± \sqrt{σ^{2} + 4σ}}{2} .

It is convenient to write

σ + 2 = 2 cos α so that p_{±} = cos α ± i sin α = e^{±iα} .

Then the boundary conditions X_{0} = 0 = X_{N} require that

A + B = 0,  and A(p_{+})^{N} + B(p_{−})^{N} = 0  which imply e^{2iNα} = 1.

So

α = \frac{mπ}{N} = mπh  for m = 1, 2 …  and σ = 2(cos α − 1) = −4 sin^{2} (\frac{mπh}{2}) = σ_{m}

say. Then X_{n} is given by

X_{n} = A (e^{inmπh} − e^{−inmπh}) = 2Ai sin(nmπh) .

For each such permissible value σ = σ_{m}, T^{j} can be found from (4.6)

T^{j+1} = λ_{m} T^{j} , where λ_{m} = [\frac{1 + σ_{m}r(1 − θ)}{1 − σ_{m}rθ}] ,

so that

T^{j} = C(λ_{m})^{j} .

Once again, we obtain the general solution by taking an arbitrary linear sum, to find

U_{n}^{j} = ∑_{m=1}^{∞} B_{m} sin(nmπh)(λ_{m})^{j} ,  (4.7)

where the coefficients B_{m} are again found from the Fourier expansion of u_{0}(x) in (4.5). We can now compare (4.7) with the exact solution (4.4) evaluated at the grid points

u_{n}^{j} ≡ u(nh, jk) = ∑_{m=1}^{∞} B_{m} sin(nmπh)(e^{−m^{2}π^{2}k})^{j}  (4.8)

The accuracy of the F.D. approximation may thus be determined by comparing λ_{m} and exp(−m^{2}π^{2}k).

Now writing β = mπh/2, we have

λ_{m} = \frac{1 − 4r(1 − θ) sin^{2} β}{1 + 4rθ sin^{2} β} = 1 − \frac{4r sin^{2} β}{1 + 4rθ sin^{2} β} .

The stability requirement is that |λ_{m}| ⩽ 1 for all m. Clearly λ_{m} ⩽ 1, while λ_{m} ⩾ −1 if

1 + 4rθ sin^{2} β ⩾ 2r sin^{2} β ,

or

(1 − 2θ)2r sin^{2} β ⩽ 1 .  (4.9)

If θ ⩾ \frac{1}{2}, therefore, the FDM (4.1) is unconditionally stable. If on the other hand 0 ⩽ θ < \frac{1}{2}, stability for all m and h can be guaranteed only if

r ⩽ \frac{1}{2(1 − 2θ)} ,  (4.10)

and (4.3) is conditionally stable. Note that when θ = 0 we recover the relation r ⩽ \frac{1}{2}.

We see that it is the high modes, the large values of m for which m ∼ N and β ∼ \frac{1}{2}π which are most likely to be unstable. This is typical behaviour for FDMs; instabilities tend to manifest themselves on the scale of the grid. The low modes, for which m = O(1), and β is small are modelled well. For them we may approximate sin β ≈ β − \frac{1}{6}β^{3} to give

λ_{m} ≈ 1 − \frac{4r(β − \frac{1}{6}β^{3} + O(β^{5}))^{2}}{1 + 4rθβ^{2} + O(rβ^{4})}

= 1 − 4rβ^{2}(1 − \frac{1}{3}β^{2})(1 − 4rθβ^{2}) + O(r^{3}β^{6})

= 1 − m^{2}π^{2}k + (θ + \frac{1}{12r})m^{4}π^{4}k^{2} + O(m^{6}k^{3})

while e^{−m^{2}π^{2}k} = 1 − m^{2}π^{2}k + \frac{1}{2}m^{4}π^{4}k^{2} + O(m^{6}k^{3}) .

The agreement between λ_{m} and exp(−m^{2}π^{2}k) is quite good, while the greatest accuracy is achieved if

θ = \frac{1}{2} − \frac{1}{12r} or r(1 − 2θ) = \frac{1}{6} .

We see from (4.10) that such a scheme would be stable. It is the generalisation of Milne’s method (θ = 0, r = \frac{1}{6}). If r is large, then the Crank-Nicolson scheme (θ = \frac{1}{2}) is close to optimal.

Most importantly, using implicit methods we can produce stable schemes with k ∼ h, and hence large values of r = k/h^{2}. Compared with the explicit schemes, we may choose relatively large time-steps.