Cost -v- benefit

7 Cost -v- benefit

At a first reading of this Section, it might be tempting to think that the extra effort involved in using Crank-Nicolson (we have to store a set of simultaneous equations, we have to solve them and we have to do this at every time-step) is enough to make the explicit method the winner in a cost-benefit analysis. But this would be wrong.

In practical problems involving numerical approximations to parabolic problems the explicit method is rarely good enough. The stability constraint ( $r \leq \frac{1}{2}$ ) imposes such tiny time-steps that it takes a great deal of time for a computer to produce approximations corresponding to even fairly modest values of $t$ . If efficiency is what matters, then Crank-Nicolson beats the explicit approach, and it is worth the extra initial effort formulating a solver (such as those we saw in HELM booklet 30) for the system of equations.

Exercises

Consider the function u defined by
$u (x, t) = x^{3} cos (x t)$

Using increments of $δ x = 0.005$ and $δ t = 0.01$ , and working to 8 decimal places, approximate
1. $u_{x} (2, 3)$ with a one-sided forward difference
2. $u_{x x} (2, 3)$ with a central difference
3. $u_{t} (2, 3)$ with a one-sided forward difference
4. $u_{t} (2, 3)$ with a central difference.
State the approximate derivatives to 3 decimal places.
The temperature $u (x, t)$ of a metal bar of length $&ell; = 3$ at a distance $x$ from one end and at time $t$ is modelled by the partial differential equation

$u_{t} = α u_{x x} (0 < x < &ell;, t > 0)$

It is given that the metal has diffusivity $α = 1.6$ , that the two ends of the bar are kept at temperature $u = 0$ and that the initial temperature distribution is

$u (x, 0) = f (x) = x (&ell; - x)$

Use the explicit difference scheme with $δ x = 0.75$ and $δ t = 0.08$ to approximate $u (x, t)$ at $t = δ t$ and $t = 2 δ t$ .
The temperature $u (x, t)$ of a metal bar of length $&ell; = 1.2$ at a distance $x$ from one end and at time $t$ is modelled by the partial differential equation
$u_{t} = α u_{x x} (0 < x < &ell;, t > 0) .$

It is given that the metal has diffusivity $α = 2.25$ , that the two ends of the bar are kept at temperature $u = 0$ and that the initial temperature distribution is

$u (x, 0) = f (x) = sin (π x ∕ &ell;)$

Use the Crank-Nicolson difference scheme with $δ x = 0.4$ and $δ t = 0.06$ to approximate $u (x, t)$ at $t = δ t$ and at $t = 2 δ t$ .

The evaluations of u we will need are u ( x , t ) = − 0.41614684 , u ( x + δ x , t ) = − 0.43162908 , u ( x − δ x , t ) = − 0.40095819 , u ( x , t + δ t ) = − 0.42521885 , u ( x , t − δ t ) = − 0 . 40703321 . It follows that
1. $u_{x} (1, 2) \approx \frac{- 0.43162908 + 0.41614684}{0.005} = - 3.096$
2. $u_{x x} (1, 2) \approx \frac{- 0.43162908 + 2 \times 0.41614684 - 0.40095819}{0.00 5^{2}} = - 11.744$
3. $u_{t} (1, 2) \approx \frac{- 0.42521885 + 0.41614684}{0.01} = - 0.907$
4. $u_{t} (1, 2) \approx \frac{- 0.42521885 + 0.40703321}{2 \times 0.01} = - 0.909$
to 3 decimal places. (Workings shown to 8 decimal places.)
In this case $r = α^{2} δ t ∕ {(δ x)}^{2} = 0.227556$ so that the numerical scheme can be written
$u_{j}^{n + 1} = u_{j}^{n} + 0.227556 (u_{j - 1}^{2} - 2 u_{j}^{n} + u_{j + 1}^{n}) = 0.544889 u_{j}^{n} + 0.227556 (u_{j - 1}^{2} + u_{j + 1}^{n})$

The first stage is to use the given data to find $u_{j}^{0}$

$\begin{matrix} u_{0}^{0} & = & 0 & from the boundary condition \\ u_{1}^{0} & = & f (δ x) = f (0.75) = 1.6875 & from the initial condition \\ u_{2}^{0} & = & f (2 δ x) = f (1.5) = 2.25 & from the initial condition \\ u_{3}^{0} & = & f (3 δ x) = f (2.25) = 1.6875 & from the initial condition \\ u_{4}^{0} & = & 0 & from the boundary condition \end{matrix}$

The first timestep will find $u_{j}^{1}$ . We note that the boundary condition implies that $u_{0}^{1} = u_{4}^{1} = 0$ .

$\begin{matrix} u_{1}^{1} & = & 0.544889 u_{1}^{0} + 0.227556 (u_{0}^{0} + u_{2}^{0}) = 0.544889 \times 1.6875 + 0.227556 (0 + 2.25) = 1.4315 \\ u_{2}^{1} & = & 0.544889 u_{2}^{0} + 0.227556 (u_{1}^{0} + u_{3}^{0}) = 0.544889 \times 2.25 + 0.227556 (1.688 + 1.688) = 1.994 \\ u_{3}^{1} & = & 0.544889 u_{3}^{0} + 0.227556 (u_{2}^{0} + u_{4}^{0}) = 0.544889 \times 1.6875 + 0.227556 (2.25 + 0) = 1.4315 \end{matrix}$

The second timestep will find $u_{j}^{2}$ . First we note that the boundary condition implies that $u_{0}^{2} = u_{4}^{2} = 0$ .

$\begin{matrix} u_{1}^{2} & = & 0.544889 u_{1}^{1} + 0.227556 (u_{0}^{1} + u_{2}^{1}) = 0.544889 \times 1.4315 + 0.227556 (0 + 1.994) = 1.233754 \\ u_{2}^{2} & = & 0.544889 u_{2}^{1} + 0.227556 (u_{1}^{1} + u_{3}^{1}) = 0.544889 \times 1.994 + 0.227556 (1.432 + 1.432) = 1.738 \\ u_{3}^{2} & = & 0.544889 u_{3}^{1} + 0.227556 (u_{2}^{1} + u_{4}^{1}) = 0.544889 \times 1.4315 + 0.227556 (1.994 + 0) = 1.233754 \end{matrix}$

where some quantities have been rounded to 6 decimal places.
In this case $r = α δ t ∕ {(δ x)}^{2} = 0.84375$ so that the numerical scheme can be written
$u_{j}^{n + 1} = u_{j}^{n} + \frac{0.84375}{2} (u_{j - 1}^{n} - 2 u_{j}^{n} + u_{j + 1}^{n} + u_{j - 1}^{n + 1} - 2 u_{j}^{n + 1} + u_{j + 1}^{n + 1})$

Moving the unknowns to the left of the equation we obtain

$- 0.42188 u_{j - 1}^{n + 1} + 1.84375 u_{j}^{n + 1} - 0.42188 u_{j + 1}^{n + 1} = 0.15625 u_{j}^{n} + 0.42188 (u_{j - 1}^{n} + u_{j + 1}^{n})$

The first stage is to use the given data to find $u_{j}^{0}$

$\begin{matrix} u_{0}^{0} & = & 0 & from the boundary condition \\ u_{1}^{0} & = & f (δ x) = f (0.4) = 0.86603 & from the initial condition \\ u_{2}^{0} & = & f (2 δ x) = f (0.8) = 0.86603 & from the initial condition \\ u_{3}^{0} & = & 0 & from the boundary condition \end{matrix}$

The first time-step will find $u_{j}^{1}$ . First we note that the boundary condition implies that $u_{0}^{1} = u_{3}^{1} = 0$ . Two uses of the stencil give

$\begin{matrix} - 0.42188 u_{0}^{1} + 1.84375 u_{1}^{1} - 0.42188 u_{2}^{1} = 0.15625 u_{1}^{0} + 0.42188 (u_{0}^{0} + u_{2}^{0}) & = & 0.50067 \\ - 0.42188 u_{1}^{1} + 1.84375 u_{2}^{1} - 0.42188 u_{3}^{1} = 0.15625 u_{2}^{0} + 0.42188 (u_{1}^{0} + u_{3}^{0}) & = & 0.50067 \end{matrix}$

The implicit nature of this method means that we have to do some extra work to complete the time-step. We must now solve the simultaneous equations

$(\begin{matrix} 1.84375 & - 0.42188 \\ - 0.42188 & 1.84375 \end{matrix}) (\begin{matrix} u_{1}^{1} \\ u_{2}^{1} \end{matrix}) = (\begin{matrix} 0.50067 \\ 0.50067 \end{matrix})$

In this case there are only two unknowns and it is a simple matter to solve the pair of equations to give $u_{1}^{1} = 0.35212$ and $u_{2}^{1} = 0.35212$ .

The second time-step will find $u_{j}^{2}$ . First we note that the boundary condition implies that $u_{0}^{2} = u_{3}^{2} = 0$ . Two uses of the stencil give

$\begin{matrix} - 0.42188 u_{0}^{2} + 1.84375 u_{1}^{2} - 0.42188 u_{2}^{2} = 0.15625 u_{1}^{1} + 0.42188 (u_{0}^{1} + u_{2}^{1}) & = & 0.20357 \\ - 0.42188 u_{1}^{2} + 1.84375 u_{2}^{2} - 0.42188 u_{3}^{2} = 0.15625 u_{2}^{1} + 0.42188 (u_{1}^{1} + u_{3}^{1}) & = & 0.20357 \end{matrix}$

The implicit nature of this method means that we have to do some extra work to complete the time-step. We must now solve the simultaneous equations

$(\begin{matrix} 1.84375 & - 0.42188 \\ - 0.42188 & 1.84375 \end{matrix}) (\begin{matrix} u_{1}^{2} \\ u_{2}^{2} \end{matrix}) = (\begin{matrix} 0.20357 \\ 0.20357 \end{matrix})$

In this case there are only two unknowns and it is a simple matter to solve the pair of equations to give $u_{1}^{2} = 0.14317$ and $u_{2}^{2} = 0.14317$ .

7 Cost -v- benefit

Exercises

Answer