The method of lines Tobin A. Driscoll Richard J. Braun 
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 Example 11.1.3 
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 x x u_t = u_{xx} u t  = u xx  periodic end conditions . Specifically, we will solve the PDE  over 0 ≤ x < 1 0\le x < 1 0 ≤ x < 1 
u ( x + 1 , t ) = u ( x , t ) for all  x . u(x+1,t)=u(x,t) \quad \text{for all $x$}. u ( x + 1 , t ) = u ( x , t ) for all  x . This is a little different from simply u ( 1 , t ) = u ( 0 , t ) u(1,t)=u(0,t) u ( 1 , t ) = u ( 0 , t ) Figure 11.2.1 
Figure 11.2.1: Left: A function whose values are the same at the endpoints of an interval does not necessarily extend to a smooth periodic function. Right: For a truly periodic function, the function values and all derivatives match at the endpoints of one period.
11.2.1 Semidiscretization ¶ As a reminder, we use u ^ \hat{u} u ^ PDE . In order to avoid carrying along redundant information about the function, we use x i = i h x_i = ih x i  = ih i = 0 , … , m − 1 i=0,\ldots,m-1 i = 0 , … , m − 1 h = 1 / m h=1/m h = 1/ m x m x_m x m  x 0 x_0 x 0  
u ^ ( x i , t ) = u ^ ( x ( i   m o d   m ) , t )   \hat{u}(x_i,t) = \hat{u}\bigl(x_{(i \bmod{m})},t \bigr) u ^ ( x i  , t ) = u ^ ( x ( i mod m )  , t ) for the exact solution u ^ \hat{u} u ^ i i i 
Next we define a vector u \mathbf{u} u 
u ( t ) = [ u 0 ( t ) u 1 ( t ) ⋮ u n ( t ) ] . \mathbf{u}(t) = \begin{bmatrix} u_0(t) \\ u_1(t) \\ \vdots \\ u_n(t) \end{bmatrix}. u ( t ) = ⎣ ⎡  u 0  ( t ) u 1  ( t ) ⋮ u n  ( t )  ⎦ ⎤  . This step is called semidiscretization , since space is discretized but time is not. As in Chapter 10 , we will replace u x x u_{xx} u xx  u \mathbf{u} u D x x \mathbf{D}_{xx} D xx  (5.4.9) (11.2.2) 
D x x = 1 h 2 [ − 2 1 1 1 − 2 1 ⋱ ⋱ ⋱ 1 − 2 1 1 1 − 2 ] . \mathbf{D}_{xx} = \frac{1}{h^2}
  \begin{bmatrix}
  -2 & 1 & & & 1 \\
  1 & -2 & 1 & & \\
  & \ddots & \ddots & \ddots & \\
  & & 1 & -2 & 1 \\
  1 &  & & 1 & -2
  \end{bmatrix}. D xx  = h 2 1  ⎣ ⎡  − 2 1 1  1 − 2 ⋱  1 ⋱ 1  ⋱ − 2 1  1 1 − 2  ⎦ ⎤  . 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 D x \mathbf{D}_x D x  
1
 2
 3
 4
 5
 6
 7
 8
 9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 """
    diffper(n, xspan)
Construct 2nd-order differentiation matrices for functions with
periodic end conditions, using `n` unique nodes in the interval
`xspan`. Returns a vector of nodes and the matrices for the first
and second derivatives.
"""
function diffper(n, xspan)
    a, b = xspan
    h = (b - a) / n
    x = @. a + h * (0:n-1)   # nodes, omitting the repeated data
    # Construct Dx by diagonals, then correct the corners.
    dp = fill(0.5 / h, n-1)       # superdiagonal
    dm = fill(-0.5 / h, n-1)      # subdiagonal
    Dx = diagm(-1 => dm, 1 => dp)
    Dx[1, n] = -1 / 2h
    Dx[n, 1] = 1 / 2h
    # Construct Dxx by diagonals, then correct the corners.
    d0 = fill(-2 / h^2, n)        # main diagonal
    dp = ones(n-1) / h^2          # superdiagonal and subdiagonal
    Dxx = diagm(-1 => dp, 0 => d0, 1 => dp)
    Dxx[1, n] = 1 / h^2
    Dxx[n, 1] = 1 / h^2
    return x, Dx, Dxx
end1
 2
 3
 4
 5
 6
 7
 8
 9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 function [x, Dx, Dxx] = diffper(n, xspan)
    %DIFFPER   Differentiation matrices for periodic end conditions. 
    % Input:
    %   n      number of subintervals (integer)
    %   xspan  endpoints of domain (vector)
    % Output:
    %   x    equispaced nodes (length n)
    %   Dx   matrix for first derivative (n by n)
    %   Dxx  matrix for second derivative (n by n)
    a = xspan(1);  b = xspan(2);
    h = (b - a) / n;
    x = a + h * (0:n-1)';   % nodes, omitting the repeated data
    % Construct Dx by diagonals, then correct the corners.
    dp =  0.5 * ones(n-1, 1) / h;    % superdiagonal
    dm = -0.5 * ones(n-1, 1) / h;    % subdiagonal
    Dx = diag(dm, -1) + diag(dp, 1);
    Dx(1, n) = -1 / (2*h);
    Dx(n, 1) =  1 / (2*h);
    % Construct Dxx by diagonals, then correct the corners.
    d0 =  -2 * ones(n, 1) / h^2;    % main diagonal
    dp =  ones(n-1, 1) / h^2;       % superdiagonal
    dm =  dp;                       % subdiagonal
    Dxx = diag(dm, -1) + diag(d0) + diag(dp, 1);
    Dxx(1, n) = 1 / h^2;
    Dxx(n, 1) = 1 / h^2;
end1
 2
 3
 4
 5
 6
 7
 8
 9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 def diffper(n, xspan):
    """
    diffper(n, xspan)
    Construct 2nd-order differentiation matrices for functions with periodic end
    conditions, using `n` unique nodes in the interval `xspan`. Return a vector of
    nodes and the  matrices for the first and second derivatives.
    """
    a, b = xspan
    h = (b - a) / n
    x = a + h * np.arange(n)  # nodes, omitting the repeated data
    # Construct Dx by diagonals, then correct the corners.
    dp = 0.5 / h * np.ones(n - 1)  # superdiagonal
    dm = -0.5 / h * np.ones(n - 1)  # subdiagonal
    Dx = np.diag(dm, -1) + np.diag(dp, 1)
    Dx[0, -1] = -1 / (2 * h)
    Dx[-1, 0] = 1 / (2 * h)
    # Construct Dxx by diagonals, then correct the corners.
    d0 = -2 / h**2 * np.ones(n)  # main diagonal
    dp = np.ones(n - 1) / h**2  # superdiagonal and subdiagonal
    Dxx = np.diag(d0) + np.diag(dp, -1) + np.diag(dp, 1)
    Dxx[0, -1] = 1 / (h**2)
    Dxx[-1, 0] = 1 / (h**2)
    return x, Dx, DxxThe PDE  u t = u x x u_t=u_{xx} u t  = u xx  
d u ( t ) d t = D x x u ( t ) ,   \frac{d \mathbf{u}(t)}{d t} = \mathbf{D}_{xx} \mathbf{u}(t), d t d u ( t )  = D xx  u ( t ) , which is simply a linear, constant-coefficient system of ordinary  differential equations. Given the initial values u ( 0 ) \mathbf{u}(0) u ( 0 ) u ( x i , 0 ) u(x_i,0) u ( x i  , 0 ) 
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.
Suppose we solve (11.2.5) IVP  integrator (6.2.5) Euler’s method  (and also AB1 from Multistep methods ). We select a time step τ \tau τ t j = j τ t_j=j\tau t j  = j τ j = 0 , 1 , … , n j=0,1,\ldots,n j = 0 , 1 , … , n u \mathbf{u} u u j ≈ u ( t j ) \mathbf{u}_j \approx \mathbf{u}(t_j) u j  ≈ u ( t j  ) j j j u j \mathbf{u}_j u j  u j u_j u j  
Thus, a fully discrete method for the heat equation is
u j + 1 = u j + τ ( D x x u j ) = ( I + τ D x x ) u j . \mathbf{u}_{j+1} = \mathbf{u}_j + \tau ( \mathbf{D}_{xx} \mathbf{u}_j) = (\mathbf{I} + \tau \mathbf{D}_{xx} ) \mathbf{u}_j. u j + 1  = u j  + τ ( D xx  u j  ) = ( I + τ D xx  ) u j  . Example 11.2.2  (Forward Euler for the heat equation)
Example 11.2.2 Let’s implement the method of Example 11.2.1 
m = 100
x, Dx, Dxx = FNC.diffper(m, [0, 1]);
tfinal = 0.15 
n = 2400           # number of time steps
τ = tfinal / n     # time step    
t = τ * (0:n)      # time valuesNext we set an initial condition. It isn’t mathematically periodic, but the end values and derivatives are so small that for numerical purposes it may as well be.
using Plots
U = zeros(m, n+1);
U[:, 1] = @. exp( -60 * (x - 0.5)^2 )
plot(x, U[:, 1];
    xaxis=(L"x"),  yaxis=(L"u(x,0)"),
    title="Initial condition")The Euler time stepping simply multiplies u j \mathbf{u}_j u j  (11.2.6) m m m 
using SparseArrays
A = sparse(I + τ * Dxx)
for j in 1:n
    U[:, j+1] = A * U[:, j]
end
plot_idx = 1:10:31
plot_times = round.(t[plot_idx], digits=4)
labels = ["t = $t" for t in plot_times]
plot(x, U[:, plot_idx];
    label=reshape(labels, 1, :),  legend=:topleft,  
    title="Heat equation by forward Euler",
    xaxis=(L"x"),  yaxis=(L"u(x,0)", [-0.25, 1]))Things seem to start well, with the initial peak widening and shrinking. But then there is a nonphysical growth in the solution.
anim = @animate for j in 1:101
    plot(x, U[:, j];
    label=@sprintf("t=%.5f", t[j]),
    xaxis=(L"x"),  yaxis=(L"u(x,t)", [-1, 2]),
    dpi=150,  title="Heat equation by forward Euler")
end
mp4(anim, "figures/diffusionFE.mp4")The growth in norm is exponential in time.
M = vec( maximum(abs, U, dims=1) )   
plot(t[1:1000], M[1:1000];
    xaxis=(L"t"),  yaxis=(:log10, L"\max_x |u(x,t)|"),
    title="Nonphysical growth")Example 11.2.2 Let’s implement the method of Example 11.2.1 
m = 100;  
[x, Dx, Dxx] = diffper(m, [0, 1]);
Ix = eye(m);
Next we set an initial condition. It isn’t mathematically periodic, but the end values and derivatives are so small that for numerical purposes it may as well be.
tfinal = 0.15;  n = 2400;  
tau = tfinal / n;  t = tau * (0:n)';
U = zeros(m, n+1);
U(:, 1) = exp( -60*(x - 0.5).^2 );
The Euler time stepping simply multiplies the solution vector by the constant matrix in (11.2.6) m m m 
A = sparse(Ix + tau * Dxx);
for j = 1:n
    U(:, j+1) = A * U(:,j);
end
index_times = 1:10:31;
show_times = t(index_times);
clf
for j = index_times
    str = sprintf("t = %.2e", t(j));
    plot(x, U(:, j), displayname=str) 
    hold on
end
legend(location="northwest")
xlabel('x'), ylabel('u(x,t)')
title('Heat equation by forward Euler')You see above that things seem to start well, with the initial peak widening and shrinking. But then there is a nonphysical growth in the solution.
clf
index_times = 1:101;
plot(x, U(:, 1))
hold on,  grid on
axis([0, 1, -1, 2])
title('Heat equation by forward Euler') 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/diffusionFE.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for frame = index_times
    cla, plot(x, U(:, frame))
    str = sprintf("t = %.3f", t(frame));
    text(0.05, 0.92, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)The growth in norm is exponential in time.
M = max(abs(U), [], 1);     % max in each column
clf,  semilogy(t, M)
xlabel('t'), ylabel('max_x |u(x,t)|') 
title('Nonphysical growth')Example 11.2.2 Let’s implement the method of Example 11.2.1 
m = 100
x, Dx, Dxx = FNC.diffper(m, [0, 1])
tfinal = 0.15  
n = 2400                 # number of time steps  
tau = tfinal / n         # time step
t = tau * arange(n+1)    # time values
Next we set an initial condition. It isn’t mathematically periodic, but the end values and derivatives are so small that for numerical purposes it may as well be.
U = zeros([m, n+1])
U[:, 0] = exp(-60 * (x - 0.5) ** 2)
plot(x, U[:, 0])
xlabel("x");  ylabel("u(x,0)")
title("Initial condition");The Euler time stepping simply multiplies the solution vector by the constant matrix in (11.2.6) m m m 
import scipy.sparse as sp
I = sp.eye(m)
A = I + tau * sp.csr_array(Dxx)
for j in range(n):
    U[:, j+1] = A @ U[:, j]
plot(x, U[:, :31:10])
ylim([-0.25, 1])
xlabel("$x$");  ylabel("$u(x,t)$")
legend([f"$t={tj:.2e}$" for tj in t[:60:20]])
title("Heat equation by forward Euler");You see above that things seem to start well, with the initial peak widening and shrinking. But then there is a nonphysical growth in the solution.
from matplotlib import animation
fig = figure()
ax = fig.add_subplot(autoscale_on=False, xlim=(0, 1), ylim=(-1, 2))
ax.grid()
line, = ax.plot([], [], '-', lw=2)
ax.set_title("Heat equation by forward Euler")
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
def animate(j):
    line.set_data(x, U[:, j])
    time_text.set_text(f"t = {t[j]:.2e}")
    return line, time_text
anim = animation.FuncAnimation(
    fig, animate, frames=range(0, 100), blit=True)
anim.save("figures/diffusionFE.mp4", fps=24)
close()
The growth in norm is exponential in time.
M = abs(U).max(axis=0)  # max in each column
semilogy(t[:500], M[:500])
xlabel("$t$");  ylabel("$\\max_x |u(x,t)|$")
title("Nonphysical growth");The method in Example 11.2.1 Example 11.2.2 Black–Scholes equation . By changing the time integrator, we can get much better results.
An alternative time discretization of (11.2.5) 
u j + 1 = u j + τ ( D x x u j + 1 ) ( I − τ D x x ) u j + 1 = u j . \begin{split}
    \mathbf{u}_{j+1} &= \mathbf{u}_j + \tau (\mathbf{D}_{xx} \mathbf{u}_{j+1})\\
    (\mathbf{I} - \tau \mathbf{D}_{xx}) \mathbf{u}_{j+1} &= \mathbf{u}_j.
\end{split} u j + 1  ( I − τ D xx  ) u j + 1   = u j  + τ ( D xx  u j + 1  ) = u j  .  Because backward Euler is an implicit method, a linear system must be solved for u j + 1 \mathbf{u}_{j+1} u j + 1  
Example 11.2.4  (Backward Euler for the heat equation)
Example 11.2.4 Now we apply backward Euler to the heat equation. We will reuse the setup from Example 11.2.2 (11.2.7) 
using SparseArrays
B = sparse(I - τ * Dxx)
factor = lu(B)
for j in 1:n
    U[:, j+1] = factor \ U[:, j]
end
using Plots
idx = 1:600:n+1
times = round.(t[idx], digits=4)
label = reshape(["t = $t" for t in times], 1, length(idx))
plot(x,U[:, idx];
    label, legend=:topleft,
    title="Heat equation by backward Euler",
    xaxis=(L"x"),  yaxis=(L"u(x,0)", [0, 1]))anim = @animate for j in 1:20:n+1
    plot(x, U[:, j];
    label=@sprintf("t=%.5f", t[j]),
    xaxis=(L"x"),  yaxis=(L"u(x,t)", [0, 1]),
    dpi=150,  title="Heat equation by backward Euler")
end
mp4(anim, "figures/diffusionBE.mp4")This solution looks physically plausible, as the large concentration in the center diffuses outward until the solution is essentially constant. Observe that the solution remains periodic in space for all time.
Example 11.2.4 Now we apply backward Euler to the heat equation. Mathematically this means multiplying by the inverse  of a matrix, but we interpret that numerically as a linear system solution. We will reuse the setup from Example 11.2.2 
B = sparse(Ix - tau * Dxx);
[l, u] = lu(B);
for j = 1:n
    U(:, j+1) = u \ (l \ U(:, j));
end
index_times = 1:600:n+1;
show_times = t(index_times);
clf
for j = index_times
    str = sprintf("t = %.2e", t(j));
    plot(x, U(:, j), displayname=str) 
    hold on
end
legend(location="northwest")
xlabel('x'), ylabel('u(x,t)')
title('Heat equation by backward Euler')clf
index_times = 1:24:n+1;
plot(x, U(:, 1))
hold on,  grid on
axis([0, 1, -0.25, 1])
title('Heat equation by backward Euler') 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/diffusionBE.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for frame = index_times
    cla, plot(x, U(:, frame))
    str = sprintf("t = %.3f", t(frame));
    text(0.05, 0.92, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)This solution looks physically plausible, as the large concentration in the center diffuses outward until the solution is essentially constant. Observe that the solution remains periodic in space for all time.
Example 11.2.4 Now we apply backward Euler to the heat equation. Mathematically this means multiplying by the inverse  of a matrix, but we interpret that numerically as a linear system solution. We will reuse the setup from Example 11.2.2 
from scipy.sparse.linalg import spsolve
B = sp.csr_matrix(I - tau * Dxx)
for j in range(n):
    U[:, j + 1] = spsolve(B, U[:, j])
plot(x, U[:, ::500])
xlabel("$x$")
ylabel("$u(x,t)$")
legend([f"$t={tj:.2g}$" for tj in t[::500]])
title("Heat equation by backward Euler");fig = figure()
ax = fig.add_subplot(autoscale_on=False, xlim=(0, 1), ylim=(-0.25, 1))
ax.grid()
line, = ax.plot([], [], '-', lw=2)
ax.set_title("Backward Euler")
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
anim = animation.FuncAnimation(
    fig, animate, frames=range(0, n+1, 20), blit=True)
anim.save("figures/diffusionBE.mp4", fps=30)
close()
This solution looks physically plausible, as the large concentration in the center diffuses outward until the solution is essentially constant. Observe that the solution remains periodic in space for all time.
Example 11.2.4 
11.2.2 Black-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) 
Example 11.2.5  (Adaptive time stepping for the heat equation)
Example 11.2.5 We set up the semidiscretization and initial condition in x x x 
m = 100
x, Dx, Dxx = FNC.diffper(m, [0, 1])
u0 = @. exp( -60*(x - 0.5)^2 );
Now, however, we apply Function 6.5.2 rk23) to the initial-value problem u ′ = D x x u \mathbf{u}'=\mathbf{D}_{xx}\mathbf{u} u ′ = D xx  u 
using OrdinaryDiffEq
tfinal = 0.25
ODE = (u, p, t) -> Dxx * u  
IVP = ODEProblem(ODE, u0, (0, tfinal))
t, u = FNC.rk23(IVP, 1e-5);
We check that the resulting solution looks realistic.
plt = plot(
    title="Heat equation by rk23",
    legend=:topleft,  
    xaxis=(L"x"),  yaxis=(L"u(x,0)", [0, 1]))
for idx in 1:600:n+1
    plot!(x, u[idx]; label="t = $(round.(t[idx], digits=4))")
end
pltanim = @animate for j in 1:20:1600
    plot(x, u[j];
    label=@sprintf("t=%.4f", t[j]),
      xaxis=(L"x"),  yaxis=(L"u(x,t)", [0, 1]),
      dpi=150,  title="Heat equation by rk23")
end
mp4(anim, "figures/diffusionRK23.mp4")The solution appears to be correct. But the number of time steps that were selected automatically is surprisingly large, considering how smoothly the solution changes.
println("Number of time steps for rk23: $(length(t)-1)")Number of time steps for rk23: 3975
 Now we apply a solver from DifferentialEquations.
u = solve(IVP, Rodas4P());
println("Number of time steps for Rodas4P: $(length(u.t) - 1)")Number of time steps for Rodas4P: 23
 The number of steps selected is reduced by a factor of more than 100!
Example 11.2.5 We set up the semidiscretization and initial condition in x x x 
m = 100;  
[x, Dx, Dxx] = diffper(m, [0, 1]);
Ix = eye(m);
u0 = exp( -60 * (x - 0.5).^2 );
Now, however, we apply a standard solver called ode45 to the initial-value problem u ′ = D x x u \mathbf{u}'=\mathbf{D}_{xx}\mathbf{u} u ′ = D xx  u 
tfinal = 0.05;
f = @(t, u) Dxx * u;
sol = ode45(f, [0, tfinal], u0);
u = @(t) deval(sol, t);
clf
for t = linspace(0, 0.05, 5)
    str = sprintf("t = %.3f", t);
    plot(x, u(t), displayname=str)
    hold on
end
xlabel("x"),  ylabel("u(x,t)")
legend()
title("Heat equation by ode45")The solution appears to be correct. But the number of time steps that were selected automatically is surprisingly large, considering how smoothly the solution changes.
time_steps_ode45 = length(sol.x) - 1Now we apply a different solver called ode15s.
sol = ode15s(f, [0, tfinal], u0);
u = @(t) deval(sol, t);
time_steps_ode15s = length(sol.x) - 1The number of steps selected was reduced by a factor of 15!
Example 11.2.5 We set up the semidiscretization and initial condition in x x x 
m = 100
x, Dx, Dxx = FNC.diffper(m, [0, 1])
u0 = exp(-60 * (x - 0.5) ** 2)
Now, however, we apply a standard solver using solve_ivp to the initial-value problem u ′ = D x x u \mathbf{u}'=\mathbf{D}_{xx}\mathbf{u} u ′ = D xx  u 
from scipy.integrate import solve_ivp
tfinal = 0.05
f = lambda t, u: Dxx @ u
sol = solve_ivp(f, [0, tfinal], u0, method="RK45", dense_output=True)
t = linspace(0, 0.05, 5)
plot(x, sol.sol(t))
xlabel("$x$"),  ylabel("$u(x,t)$")
legend([f"$t={tj:.4g}$" for tj in t])
title("Heat equation by RK45");The solution appears to be correct. But the number of time steps that were selected automatically is surprisingly large, considering how smoothly the solution changes.
print(f"RK45 took {len(sol.t) - 1} steps")Now we apply a different solver called BDF.
sol = solve_ivp(f, [0, tfinal], u0, method="BDF")
print(f"BDF took {len(sol.t) - 1} steps")The number of steps selected was reduced by a factor of 20!
The adaptive time integrators can all produce solutions. But, as seen in Example 11.2.5 (11.2.5) 
11.2.3 Exercises ¶ ⌨ Revisit Example 11.2.2 m = 20 , 30 , … , 120 m=20,30,\dots,120 m = 20 , 30 , … , 120 n = 20 , 30 , 40 , … n=20,30,40,\dots n = 20 , 30 , 40 , … n n n N ( m ) N(m) N ( m ) N ( m ) N(m) N ( m ) m m m N = O ( m p ) N=O(m^p) N = O ( m p ) p p p p p p 
In Example 11.2.5 t → ∞ t\to \infty t → ∞ u ( x , t ) u(x,t) u ( x , t ) 
(a)  ⌨ Set m = 400 m=400 m = 400 IVP  solver, as shown in Example 11.2.5 
(b)  ✍ Prove that Q = ∫ 0 1 u ( x , t )   d x Q = \int_0^1 u(x,t) \,dx Q = ∫ 0 1  u ( x , t ) d x t t t PDE  and periodic end conditions.)
(c)  ⌨ Use Function 5.6.1 Q Q Q t = 0 t=0 t = 0 
✍ Apply the trapezoid IVP  formula (AM2 in Table 6.6.1 (11.2.5) Crank–Nicolson  method:
( I − 1 2 τ D x x ) u j + 1 = ( I + 1 2 τ D x x ) u j . (\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}. ( I − 2 1  τ D xx  ) u j + 1  = ( I + 2 1  τ D xx  ) u j  . Note that each side of the method is evaluated at a different time level.
The PDE  u t = 2 u + u x x u_t = 2u + u_{xx} u t  = 2 u + u xx  
(a)  ✍ Derive an equation analogous to (11.2.7) 
(b)  ⌨ Apply your formula from part (a) to solve this PDE  with periodic boundary conditions for the same initial condition as in Example 11.2.4 m = 200 m=200 m = 200 n = 1000 n=1000 n = 1000 t = 0 , 0.04 , 0.08 , … , 0.2 t=0,0.04,0.08,\ldots,0.2 t = 0 , 0.04 , 0.08 , … , 0.2 0 ≤ t ≤ 0.2 0\le t \le 0.2 0 ≤ t ≤ 0.2 
✍ In this problem, you will analyze the convergence of the explicit method given by (11.2.6) u i , j u_{i,j} u i , j  x i x_i x i  t j t_j t j  
(a)  Write the method in scalar form as
u i , j + 1 = ( 1 − 2 λ ) u i , j + λ u i + 1 , j + λ u i − 1 , j , u_{i,j+1} = (1-2\lambda) u_{i,j} + \lambda u_{i+1,j} + \lambda u_{i-1,j}, u i , j + 1  = ( 1 − 2 λ ) u i , j  + λ u i + 1 , j  + λ u i − 1 , j  , where λ = τ / h 2 > 0 \lambda = \tau/h^2>0 λ = τ / h 2 > 0 
(b)  Taylor series of the exact solution u ^ \hat{u} u ^ 
u ^ i , j + 1 = u i , j + ∂ u ^ ∂ t ( x i , t j ) τ + O ( τ 2 ) , u ^ i ± 1 , j = u ^ i , j ± ∂ u ^ ∂ x ( x i , t j ) h + ∂ 2 u ^ ∂ x 2 ( x i , t j ) h 2 2 ± ∂ 3 u ^ ∂ x 3 ( x i , t j ) h 3 6 + O ( h 4 ) . \begin{align*}
\hat{u}_{i,j+1} &= u_{i,j} + \frac{\partial \hat{u}}{\partial t} (x_i,t_j) \tau + O(\tau^2),\\
% \frac{\partial^2 u}{\partial t^2} (x_i,\bar{t}) \frac{\tau^2}{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).
%\frac{\partial^4 u}{\partial x^4} (\bar{x}_\pm,t_j) \frac{h^4}{24}.
\end{align*} u ^ i , j + 1  u ^ i ± 1 , j   = u i , j  + ∂ t ∂ u ^  ( x i  , t j  ) τ + O ( τ 2 ) , = u ^ i , j  ± ∂ x ∂ u ^  ( x i  , t j  ) h + ∂ x 2 ∂ 2 u ^  ( x i  , t j  ) 2 h 2  ± ∂ x 3 ∂ 3 u ^  ( x i  , t j  ) 6 h 3  + O ( h 4 ) .  Use these to show that
u ^ i , j + 1 = [ ( 1 − 2 λ ) u ^ i , j + λ u ^ i + 1 , j + λ u ^ i − 1 , j ] + O ( τ 2 + h 2 ) = F ( λ , u ^ i , j , u ^ i + 1 , j , u ^ i − 1 , j ) + O ( τ 2 + h 2 ) . \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*} u ^ i , j + 1   = [ ( 1 − 2 λ ) u ^ i , j  + λ u ^ i + 1 , j  + λ u ^ i − 1 , j  ] + O ( τ 2 + h 2 ) = F ( λ , u ^ i , j  , u ^ i + 1 , j  , u ^ i − 1 , j  ) + O ( τ 2 + h 2 ) .  (The last line should be considered a definition of the function F F F 
(c)  The numerical solution satisfies
u i , j + 1 = F ( λ , u i , j , u i + 1 , j , u i − 1 , j ) u_{i,j+1}=F\bigl( \lambda,u_{i,j}, u_{i+1,j} , u_{i-1,j}\bigr) u i , j + 1  = F ( λ , u i , j  , u i + 1 , j  , u i − 1 , j  ) exactly. Using this fact, subtract u i , j + 1 u_{i,j+1} u i , j + 1  
e i , j + 1 = F ( λ , e i , j , e i + 1 , j , e i − 1 , j ) + O ( τ 2 + h 2 ) , 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), e i , j + 1  = F ( λ , e i , j  , e i + 1 , j  , e i − 1 , j  ) + O ( τ 2 + h 2 ) , where e i , j = u ^ i , j − u i , j e_{i,j}=\hat{u}_{i,j}-u_{i,j} e i , j  = u ^ i , j  − u i , j  i i i j j j 
(d)  Define E j E_j E j  ∣ e i , j ∣ |e_{i,j}| ∣ e i , j  ∣ i i i λ < 1 / 2 \lambda<1/2 λ < 1/2 h h h τ \tau τ τ \tau τ h h h 
E j + 1 = E j + O ( τ 2 + h 2 ) ≤ E j + K j ( τ 2 + h 2 ) E_{j+1}  = E_{j} + O\Bigl(\tau^2+h^2\Bigr) \le E_{j} + K_j\bigl(\tau^2+h^2\bigr) E j + 1  = E j  + O ( τ 2 + h 2 ) ≤ E j  + K j  ( τ 2 + h 2 ) for a positive K j K_j K j  τ \tau τ h h h 
(e)  If the initial conditions are exact, then E 0 = 0 E_0=0 E 0  = 0 K j K_j K j  λ < 1 / 2 \lambda<1/2 λ < 1/2 E n = O ( τ ) E_n = O(\tau) E n  = O ( τ ) τ → 0 \tau\to 0 τ → 0