The method of lines

Our strategy in Black–Scholes equation was to discretize both the time and space derivatives using finite differences, then rearrange so that we could march the solution forward through time. It was partially effective, but as Demo 11.1.3 shows, not a sure thing, for reasons we look into starting in the next section.

First, though, we want to look at a broader version of the discretization approach. To introduce ideas, let’s use the simpler heat equation, $u_t = u_{xx}$, as a model. Because boundaries always complicate things, we will start by doing the next best thing to having no boundaries at all: periodic end conditions. Specifically, we will solve the PDE over $0\le x < 1$ and require at all times that

This is a little different from simply $u(1,t)=u(0,t)$, as Figure 11.2.1 illustrates.

11.2.1Semidiscretization¶

As always, we use $\hat{u}$ when we specifically refer to the exact solution of the PDE. In order to avoid carrying along redundant information about the function, we use $x_i = ih$ only for $i=0,\ldots,m-1$, where $h=1/m$, and it’s understood that a reference to $x_m$ is silently translated to one at $x_0$. More generally, we have the identity

for the exact solution $\hat{u}$ at any value of $i$.

Next we define a vector $\mathbf{u}$ by

This step is called semidiscretization, since space is discretized but time is not. As in Chapter 10, we will replace $u_{xx}$ with multiplication of $\mathbf{u}$ by a differentiation matrix $\mathbf{D}_{xx}$. The canonical choice is the three-point finite-difference formula (5.4.12), which in light of the periodicity (11.2.2) leads to

Note well how the first and last rows have elements that “wrap around” from one end of the domain to the other by periodicity. Because we will be using this matrix quite a lot, we create Function 11.2.1 to compute it, as well as the corresponding second-order first derivative matrix $\mathbf{D}_x$ for periodic end conditions.

The PDE $u_t=u_{xx}$ is now approximated by the semidiscrete problem

which is simply a linear, constant-coefficient system of ordinary differential equations. Given the initial values $\mathbf{u}(0)$ obtained from $u(x_i,0)$, we have an initial-value problem that we already know how to solve!

Semidiscretization is often called the method of lines. Despite the name, it is not exactly a single method because both space and time discretizations have to be specified in order to get a concrete algorithm. The key concept is the separation of those two discretizations, and in that way, it’s related to separation of variables in analytic methods for the heat equation.

The method in Example 11.2.1 and Demo 11.2.2 is essentially the same one we used for the Black–Scholes equation in Black–Scholes equation. By changing the time integrator, we can get much better results.

Demo 11.2.4 suggests that implicit time stepping methods have an important role in diffusion. We will analyze the reason in the next few sections.

11.2.2Black-box IVP solvers¶

Instead of coding one of the Runge–Kutta or multistep formulas directly for a method of lines solution, we could use any of the IVP solvers from Chapter 6, or a solver from the DifferentialEquations package, to solve the ODE initial-value problem (11.2.5).

The adaptive time integrators can all produce solutions. But, as seen in Demo 11.2.5, they are not equivalent in every important sense. Whether we choose to implement a method directly with a fixed step size, or automatically with adaptation, there is something crucial to understand about the semidiscrete problem (11.2.5) that will occupy our attention in the next two sections.

11.2.3Exercises¶

1. ⌨ Revisit Demo 11.2.2. For each $m=20,30,\dots,120$ points in space, let $n=20,30,40,\dots$ in turn until you reach the smallest $n$ such that the numerical solution remains bounded above by 2 for all time; call this value $N(m)$. Make a log-log plot of $N(m)$ as a function of $m$. If you suppose that $N=O(m^p)$ for a simple rational number $p$, what is a reasonable hypothesis for $p$?

2. In Demo 11.2.5, as $t\to \infty$ the solution $u(x,t)$ approaches a value that is constant in both space and time.

   (a) ⌨ Set $m=400$ and use Rodas4P, as shown in Demo 11.2.5, to find this constant value to at least eight digits of accuracy.

   (b) ✍ Prove that $\int_0^1 u(x,t) \,dx$ is constant in time.

   (c) ⌨ Use Function 5.6.1 on the initial condition function, and compare to the result of part (a).

3. ✍ Apply the trapezoid IVP formula (AM2) to the semidiscretization (11.2.5) and derive what is known as the Crank–Nicolson method:

   $$(\mathbf{I} - \tfrac{1}{2}\tau \mathbf{D}_{xx}) \mathbf{u}_{j+1} = (\mathbf{I} + \tfrac{1}{2}\tau \mathbf{D}_{xx}) \mathbf{u}_{j}.$$

   (11.2.8)

   Note that each side of the method is evaluated at a different time level.

4. ⌨ Repeat Demo 11.2.4 using the Crank–Nicolson method (11.2.8). Then try for $n=240$ as well, which uses a time step ten times larger than before. Does the solution remain stable?

5. The PDE $u_t = 2u + u_{xx}$ combines growth with diffusion.

   (a) ✍ Derive an equation analogous to (11.2.7) that combines second-order semidiscretization in space with the backward Euler solver in time.

   (b) ⌨ Apply your formula from part (a) to solve this PDE with periodic boundary conditions for the same initial condition as in Demo 11.2.4. Use $m=200$ points in space and $n=1000$ time levels. Plot the solution on one graph at times $t=0,0.04,0.08,\ldots,0.2$, or animate the solution over $0\le t \le 0.2$.

6. ✍ In this problem, you will analyze the convergence of the explicit method given by (11.2.6). Recall that the discrete approximation $u_{i,j}$ approximates the solution at $x_i$ and $t_j$.

   (a) Write the method in scalar form as

   $$u_{i,j+1} = (1-2\lambda) u_{i,j} + \lambda u_{i+1,j} + \lambda u_{i-1,j},$$

   (11.2.9)

   where $\lambda = \tau/h^2>0$.

   (b) Taylor series of the exact solution $\hat{u}$ imply that

   $$\begin{align*} \hat{u}_{i,j+1} &= u_{i,j} + \frac{\partial \hat{u}}{\partial t} (x_i,t_j) \tau + O(\tau^2),\\ \hat{u}_{i\pm1,j} &= \hat{u}_{i,j} \pm \frac{\partial \hat{u}}{\partial x} (x_i,t_j) h + \frac{\partial^2 \hat{u}}{\partial x^2} (x_i,t_j) \frac{h^2}{2} \pm \frac{\partial^3 \hat{u}}{\partial x^3} (x_i,t_j) \frac{h^3}{6}+ O(h^4). \end{align*}$$

   (11.2.10)

   Use these to show that

   $$\begin{align*} \hat{u}_{i,j+1} & = \left[ (1-2\lambda) \hat{u}_{i,j} + \lambda \hat{u}_{i+1,j} + \lambda \hat{u}_{i-1,j}\right] + O\Bigl(\tau^2+h^2 \Bigr)\\ &= F\left( \lambda,\hat{u}_{i,j}, \hat{u}_{i+1,j} , \hat{u}_{i-1,j}\right) + O\Bigl(\tau^2+h^2\Bigr). \end{align*}$$

   (11.2.11)

   (The last line should be considered a definition of the function $F$.)

   (c) The numerical solution satisfies

   $$u_{i,j+1}=F\bigl( \lambda,u_{i,j}, u_{i+1,j} , u_{i-1,j}\bigr)$$

   (11.2.12)

   exactly. Using this fact, subtract $u_{i,j+1}$ from both sides of the last line in part (b) to show that

   $$e_{i,j+1} = F\left( \lambda,e_{i,j}, e_{i+1,j} ,e_{i-1,j}\right) + O\Bigl(\tau^2+h^2\Bigr),$$

   (11.2.13)

   where $e_{i,j}=\hat{u}_{i,j}-u_{i,j}$ is the error in the numerical solution for all $i$ and $j$ .

   (d) Define $E_j$ as the maximum of $|e_{i,j}|$ over all values of $i$, and use the result of part (c) to show that if $\lambda<1/2$ is kept fixed as $h$ and τ approach zero, then for sufficiently small τ and $h$,

   $$E_{j+1} = E_{j} + O\Bigl(\tau^2+h^2\Bigr) \le E_{j} + K_j\bigl(\tau^2+h^2\bigr)$$

   (11.2.14)

   for a positive $K_j$ independent of τ and $h$.

   (e) If the initial conditions are exact, then $E_0=0$. Use this to show finally that if the $K_j$ are bounded above and $\lambda<1/2$ is kept fixed, then $E_n = O(\tau)$ as $\tau\to 0$.