Boundaries - Fundamentals of Numerical Computation

So far we have considered the method of lines for problems with periodic end conditions, which is much like having no boundary at all. How can boundary conditions be incorporated into this technique?

Suppose we are given a nonlinear PDE of the form

u_t = \phi(t,x,u,u_x,u_{xx}), \quad a \le x \le b.

(11.5.1)

Not all such PDEs are parabolic (essentially, including diffusion), but we will assume this to be the case. Suppose also that the solution is subject to the boundary conditions

\begin{split} g_1\left( u(a,t), \frac{\partial u}{\partial x}(a,t) \right) &= 0, \\ g_2\left( u(b,t), \frac{\partial u}{\partial x}(b,t) \right) &= 0. \\ \end{split}

(11.5.2)

These include Dirichlet, Neumann, and Robin conditions, which are the linear cases of (11.5.2).

11.5.1Boundary removal¶

As usual, we replace $u(x,t)$ by the semidiscretized $\mathbf{u}(t)$ , where $u_i(t)\approx \hat{u}(x_i,t)$ and $i=0,\ldots,m$ . We require the endpoints of the interval to be included in the discretization, that is, $x_0=a$ and $x_m=b$ . Then we have a division of the semidiscrete unknown $\mathbf{u}(t)$ into interior and boundary nodes:

\mathbf{u} = \begin{bmatrix} u_0 \\ \mathbf{v} \\ u_m \end{bmatrix},

(11.5.3)

where $\mathbf{v}$ are the solution values over the interior of the interval. The guiding principle is to let the interior unknowns $\mathbf{v}$ be governed by a discrete form of the PDE, while the endpoint values are chosen to satisfy the boundary conditions. As a result, we will develop an initial-value problem for the interior unknowns only:

\frac{d \mathbf{v}}{d t} = \mathbf{f}(t,\mathbf{v}).

(11.5.4)

The boundary conditions are used only in the definition of $\mathbf{f}$ . As in Nonlinearity and boundary conditions, define

\mathbf{u}' = \mathbf{D}_x \mathbf{u}.

(11.5.5)

Then Equation (11.5.2) takes the form

\begin{split} g_1( u_0, u'_0 ) &= 0, \\ g_2( u_m, u'_m ) &= 0. \end{split}

(11.5.6)

Given a value of $\mathbf{v}$ for the interior nodes, Equation (11.5.6) may be considered a system of two equations for the unknown boundary values $u_0$ and $u_m$ . This system will be a linear one for Dirichlet, Neumann, and Robin conditions.

Example 11.5.1

Recall the Black–Scholes PDE (11.1.7),

u_t = \frac{1}{2} \sigma^2 x^2 u_{xx} + rx u_x - ru,

(11.5.7)

subject to $u(0)=0$ and $u_x(S_\text{max})=1$ . These imply $u_0=0$ and

\frac{ \tfrac{1}{2} u_{m-2} -2u_{m-1} + \tfrac{3}{2} u_{m}}{h} = 1.

(11.5.8)

Hence

\begin{bmatrix} 1 & 0 \\ 0 & \frac{3}{2} \end{bmatrix} \begin{bmatrix} u_0 \\ u_m \end{bmatrix} = \begin{bmatrix} 0 \\ h \end{bmatrix} - \begin{bmatrix} 0 & \cdots & 0 & 0 \\ 0 & \cdots & -\frac{1}{2} & 2 \end{bmatrix} \mathbf{v}.

(11.5.9)

11.5.2Implementation¶

The steps to evaluate $\mathbf{f}$ in (11.5.4) now go as follows.

Our full implementation of the method of lines for (11.5.1)--(11.5.2) is given in Function 11.5.2. It uses Function 10.3.2 (diffcheb) to set up a Chebyshev discretization. The nested function extend performs steps 1--2 of Algorithm 11.5.1 by calling Function 4.6.3 (levenberg) to solve the potentially nonlinear system (11.5.6). Then it sets up and solves an IVP, adding steps 3--4 of Algorithm 11.5.1 within the ode! function. Finally, it returns the node vector x and a function of t that applies extend to $\mathbf{v}(t)$ to compute $\mathbf{u}(t)$ .

Algorithm 11.5.2 (parabolic)

Julia

MATLAB

Python

Solution of parabolic PDEs by the method of lines

parabolic.jl

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
30
31
32
33
34
35
36
"""
    parabolic(ϕ, xspan, m, g₁, g₂, tspan, init)

Solve a parabolic PDE by the method of lines. The PDE is
∂u/∂t = `ϕ`(t, x, u, ∂u/∂x, ∂^2u/∂x^2), `xspan` gives the space
domain, m gives the degree of a Chebyshev spectral discretization,
`g₁` and `g₂` are functions of (u,∂u/∂x) at the domain ends that
should be made zero, `tspan` is the time domain, and `init` is a
function of x that gives the initial condition. Returns a vector
`x` and a function of t that gives the semidiscrete solution at `x`.
"""
function parabolic(ϕ, xspan, m, g₁, g₂, tspan, init)
    x, Dₓ, Dₓₓ = diffcheb(m, xspan)
    int = 2:m    # indexes of interior nodes

    function extend(v)
        function objective(ubc)
            u₀, uₘ = ubc
            uₓ = Dₓ * [u₀; v; uₘ]
            return [g₁(u₀, uₓ[1]), g₂(uₘ, uₓ[end])]
        end
        ubc = levenberg(objective, [0, 0])[end]
        return [ubc[1]; v; ubc[2]]
    end

    function ode!(f, v, p, t)
        u = extend(v)
        uₓ, uₓₓ = Dₓ * u, Dₓₓ * u
        @. f = ϕ(t, x[int], u[int], uₓ[int], uₓₓ[int])
    end

    ivp = ODEProblem(ode!, init.(x[int]), float.(tspan))
    u = solve(ivp)

    return x, t -> extend(u(t))
end

Solution of parabolic PDEs by the method of lines

parabolic.m

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
30
31
32
33
34
35
36
37
function [x, u] = parabolic(phi, xspan, m, ga, gb, tspan, init)
% PARABOLIC   Solve parabolic PDE by the method of lines.
% Input:
%   phi      defines ∂u/∂t = phi(t, x, u, ∂u/∂x, ∂^2u/∂x^2) 
%   xspan    spatial domain
%   m        number of spatial nodes
%   ga, gb   boundary conditions as functions of u and ∂u/∂x
%   tspan    time interval
%   init     initial condition as a function of x
% Output:       
%   x        spatial nodes (vector)
%   u        function for the solution u(t) at nodes

    [x, Dx, Dxx] = diffcheb(m, xspan);
    int = 2:m;    % indexes of interior nodes

    function u = extend(v)
        function residual = objective(ubc)
            ua = ubc(1);  ub = ubc(2);
            ux = Dx * [ua; v; ub];
            residual = [ga(ua, ux(1)); gb(ub, ux(end))];
        end
        ubc = levenberg(@objective, [0, 0]);
        ubc = ubc(:, end);
        u = [ubc(1); v; ubc(2)];
    end

    function f = mol_ode(t, v)
        u = extend(v);
        ux = Dx * u;
        uxx = Dxx * u;
        f = phi(t, x(int), u(int), ux(int), uxx(int));
    end

    sol = ode15s(@mol_ode, tspan, init(x(int)));
    u = @(t) extend(deval(sol, t));
end

Solution of parabolic PDEs by the method of lines

parabolic.py

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
def parabolic(phi, xspan, m, ga, gb, tspan, init):
    """
        parabolic(phi, xspan, m, ga, gb, tspan, init)

    Solve a parabolic PDE by the method of lines. The PDE is 
    ∂u/∂t = phi(t,x,u,∂u/∂x,∂^2u/∂x^2), xspan gives the space 
    domain, m gives the degree of a Chebyshev spectral discretization, ga and gb are functions of (u,∂u/∂x) at the domain ends that should be made zero, tspan is the time domain, and init is a function of x that gives the initial condition. Returns a vector x and a function of t that gives the semidiscrete solution at x. 
    """   
    x, Dx, Dxx = diffcheb(m, xspan)
    int = range(1, m)    # indexes of interior nodes

    def extend(v):
        def objective(ubc):
            u0, um = ubc
            ux = Dx @ np.hstack([u0, v, um])
            return np.array([ga(u0, ux[1]), gb(um, ux[-1])])
        ubc = levenberg(objective, np.array([0, 0]))[-1]
        return np.hstack([ubc[0], v, ubc[-1]])

    def ode(t, v):
        u = extend(v)
        ux, uxx = Dx @ u, Dxx @ u
        return phi(t, x[int], u[int], ux[int], uxx[int])

    u0 = init(x[int])
    sol = solve_ivp(ode, tspan, u0, method="BDF", dense_output=True)

    return x, lambda t: extend(sol.sol(t))

In many specific problems, extend does more work than is truly necessary. Dirichlet boundary conditions, for instance, define $u_0$ and $u_m$ directly, and there is no need to solve a nonlinear system. Exercise 11.5.6 asks you to modify the code to take advantage of this case. The price of solving a more general set of problems in Function 11.5.2 is some speed in such special cases.^[1]

Example 11.5.3 (Heat equation with Dirichlet boundary conditions)

Let’s solve the heat equation $u_t=u_{xx}$ on $[-1,1]$ , subject to the Dirichlet boundary conditions $u(-1,t)=0$ , $u(1,t)=2$ .

Julia

MATLAB

Python

Example 11.5.3

First, we define functions for the PDE and each boundary condition.

ϕ = (t, x, u, uₓ, uₓₓ) -> uₓₓ
g₁ = (u, uₓ) -> u
g₂ = (u, uₓ) -> u - 2;

Our next step is to write a function to define the initial condition. This one satisfies the boundary conditions exactly.

init = x -> 1 + sinpi(x/2) + 3 * (1-x^2) * exp(-4x^2);

Now we can use Function 11.5.2 to solve the problem.

using Plots
x, u = FNC.parabolic(ϕ, (-1, 1), 60, g₁, g₂, (0, 0.75), init)
plt = plot(
    xlabel=L"x",  ylabel=L"u(x,t)",
    legend=:topleft,  title="Solution of the heat equation")
for t in 0:0.1:0.4
    plot!(x, u(t), label=@sprintf("t=%.2f", t))
end
plt

Source

anim = @animate for t in range(0,0.75,length=201) 
    plot(x, u(t);
        label=@sprintf("t=%.2f", t),  legend=:topleft,
        xaxis=(L"x"),  yaxis=(L"u(x,t)", (0, 4.2)), 
        title="Heat equation",  dpi=150)
end
mp4(anim, "figures/boundaries-heat.mp4", fps=30)

Example 11.5.3

First, we define functions for the PDE and each boundary condition.

phi = @(t, x, u, ux, uxx) uxx;
ga = @(u, ux) u;
gb = @(u, ux) u - 2;

Our next step is to write a function to define the initial condition. This one satisfies the boundary conditions exactly.

init = @(x) 1 + sin(pi * x/2) + 3 * (1 - x.^2) .* exp(-4*x.^2);

Now we can use Function 11.5.2 to solve the problem.

[x, u] = parabolic(phi, [-1, 1], 60, ga, gb, [0, 0.75], init);

clf
for t = 0:0.1:0.5
    str = sprintf("t = %.2f", t);
    plot(x, u(t), displayname=str)
    hold on
end
xlabel("x"),  ylabel("u(x,t)")
legend()
title("Heat equation with Dirichlet boundaries")

clf
plot(x, u(0))
hold on,  grid on
axis([-1, 1, 0, 4.2])
title('Heat equation with Dirichlet boundaries') 
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/boundaries-heat.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for t = linspace(0, 0.75, 201)
    cla, plot(x, u(t))
    str = sprintf("t = %.3f", t);
    text(-0.9, 3.8, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

Example 11.5.3

First, we define functions for the PDE and each boundary condition.

phi = lambda t, x, u, ux, uxx: uxx
ga = lambda u, ux: u
gb = lambda u, ux: u - 2

Our next step is to write a function to define the initial condition. This one satisfies the boundary conditions exactly.

init = lambda x: 1 + sin(pi * x/2) + 3 * (1 - x**2) * exp(-4*x**2)

Now we can use Function 11.5.2 to solve the problem.

x, u = FNC.parabolic(phi, (-1, 1), 60, ga, gb, (0, 0.75), init)

for t in arange(0,0.5,0.1):
    plot(x, u(t), label=f"t={t:.2f}")
xlabel("$x$"),  ylabel("$u(x,t)$")
legend()
title("Heat equation with Dirichlet boundaries");

Source

from matplotlib import animation
fig = figure()
ax = fig.add_subplot(autoscale_on=False, xlim=(-1, 1), ylim=(0, 4.2))

line, = ax.plot([], [], '-')
ax.set_title("Heat equation with Dirichlet boundaries")
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

def animate(t):
    line.set_data(x, u(t))
    time_text.set_text(f"t = {t:.2e}")
    return line, time_text

anim = animation.FuncAnimation(
    fig, animate, frames=linspace(0, 0.75, 201), blit=True)
anim.save("figures/boundaries-heat.mp4", fps=30)
close()

Example 11.5.4 (Heat equation with nonlinear source)

We solve a heat equation with a nonlinear source term,

u_t = u^2 + u_{xx}.

(11.5.12)

One interpretation of this PDE is an exothermic chemical reaction whose rate increases with temperature. We solve over $x \in [0,1]$ with homogeneous conditions of different kinds.

Julia

MATLAB

Python

Example 11.5.4

ϕ = (t, x, u, uₓ, uₓₓ) -> u^2 + uₓₓ
g₁ = (u, uₓ) -> u
g₂ = (u, uₓ) -> uₓ
init = x -> 400x^4 * (1 - x)^2
x, u = FNC.parabolic(ϕ, (0, 1), 60, g₁, g₂, (0, 0.1), init);

Source

anim = @animate for t in range(0, 0.1, length=101) 
    plot(x, u(t);
        label=@sprintf("t=%.4f", t),  legend=:topleft,
        xaxis=(L"x"),  yaxis=(L"u(x,t)", (0, 10)),
        dpi=150, title="Heat equation with source")
end
mp4(anim, "figures/boundaries-source.mp4", fps=30)

Example 11.5.4

phi = @(t, x, u, ux, uxx) u.^2 + uxx;
ga = @(u, ux) u;
gb = @(u, ux) ux;

init = @(x) 400 * x.^4 .* (1 - x).^2;
[x, u] = parabolic(phi, [0, 1], 60, ga, gb, [0, 0.1], init);

clf
plot(x, u(0))
hold on,  grid on
axis([0, 1, 0, 10])
title("Heat equation with source")
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/boundaries-source.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for t = linspace(0, 0.1, 101)
    cla, plot(x, u(t))
    str = sprintf("t = %.3f", t);
    text(0.05, 9.2, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

Example 11.5.4

phi = lambda t, x, u, ux, uxx: u**2 + uxx
ga = lambda u, ux: u
gb = lambda u, ux: ux
init = lambda x: 400 * x**4 * (1 - x)**2
x, u = FNC.parabolic(phi, (0, 1), 60, ga, gb, (0, 0.1), init);

Source

fig = figure()
ax = fig.add_subplot(autoscale_on=False, xlim=(0, 1), ylim=(0, 10))

line, = ax.plot([], [], '-')
ax.set_title("Heat equation with source")
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

anim = animation.FuncAnimation(
    fig, animate, frames=linspace(0, 0.1, 101), blit=True)
anim.save("figures/boundaries-source.mp4", fps=30)
close()

Finally, we return to the example of the Black–Scholes equation from Black–Scholes equation.

Example 11.5.5 (Black–Scholes equation with mixed boundary conditions)

We solve the Black–Scholes PDE (11.1.2) with initial condition $u(x,0) = \max\{0,x-K\}$ and the boundary conditions $u(0,t)=0$ and $u_x(S_\text{max},t)=1$ . We choose $S_\text{max}=8$ , $r=0.08$ , $\sigma=0.06$ , and $K=3$ .

Julia

MATLAB

Python

Example 11.5.5

K = 3;  σ = 0.06;  r = 0.08;  Smax = 8;
ϕ = (t, x, u, uₓ, uₓₓ) -> σ^2/2 * (x^2 * uₓₓ) + r*x*uₓ - r*u
g₁ = (u, uₓ) -> u
g₂ = (u, uₓ) -> uₓ - 1;

u₀ = x -> max(0, x - K)
x, u = FNC.parabolic(ϕ, (0, Smax), 80, g₁, g₂, (0, 15), u₀);

Source

anim = @animate for t in range(0, 15, 151) 
    plot(x, u(t);
        label=@sprintf("t=%.4f", t),  legend=:topleft,
        xaxis=(L"x"),  yaxis=(L"u(x,t)", (-0.5, 8)), 
        dpi=150,  title="Black–Scholes equation")
end
mp4(anim, "figures/boundaries-bs.mp4", fps=30)

Recall that $u$ is the value of the call option, and time runs backward from the strike time. The longer the horizon, the more value the option has due to anticipated growth in the stock price.

Example 11.5.5

K = 3;  sigma = 0.06;  r = 0.08;  Smax = 8;
phi = @(t, x, u, ux, uxx) sigma.^2/2 * (x.^2 .* uxx) + r * x.*ux - r*u;
ga = @(u, ux) u;
gb = @(u, ux) ux - 1;

init = @(x) max(0, x - K);
[x, u] = parabolic(phi, [0, Smax], 80, ga, gb, [0, 15], init);

clf
plot(x, u(0))
hold on,  grid on
axis([0, Smax, -0.1, 8])
title("Black–Scholes equation with boundaries")
xlabel('x'),  ylabel('u(x,t)')
vid = VideoWriter("figures/boundaries-bs.mp4", "MPEG-4");
vid.Quality = 85;
open(vid);
for t = linspace(0, 15, 151)
    cla, plot(x, u(t))
    str = sprintf("t = %.1f", t);
    text(0.5, 7.1, str);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

Recall that $u$ is the value of the call option, and time runs backward from the strike time. The longer the horizon, the more value the option has due to anticipated growth in the stock price.

Example 11.5.5

K = 3;  sigma = 0.06;  r = 0.08;  Smax = 8;
phi = lambda t, x, u, ux, uxx: sigma**2/2 * (x**2 * uxx) + r*x*ux - r*u
ga = lambda u, ux: u
gb = lambda u, ux: ux - 1

u0 = lambda x: maximum(0, x - K)
x, u = FNC.parabolic(phi, (0, Smax), 80, ga, gb, (0, 15), u0);

Source

fig = figure()
ax = fig.add_subplot(autoscale_on=False, xlim=(0, Smax), ylim=(-0.5, 8))

line, = ax.plot([], [], '-')
ax.set_title("Black–Scholes equation with boundaries")
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)

anim = animation.FuncAnimation(
    fig, animate, frames=linspace(0, 15, 151), blit=True)
anim.save("figures/boundaries-bs.mp4", fps=30)
close()

Recall that $u$ is the value of the call option, and time runs backward from the strike time. The longer the horizon, the more value the option has due to anticipated growth in the stock price.

11.5.3Exercises¶

Exercise 11.5.3

Consider Example 11.5.4, combining diffusion with a nonlinear source term.

(a) ✍ Suppose we ignore the diffusion. Use separation of variables (or computer algebra) to solve the IVP $u_t=u^2$ , $u(0) = A>0$ . What happens as $t\to 1/A$ from below?

(b) ⌨ Try to continue the solution in the demo to $t=0.5$ (you may get a warning from the IVP solver). Plot the solution at $t=0.15$ , $t=0.2$ , and $t=0.25$ .

(c) ⌨ Let the initial condition be $u(x,0) = C x^4(1-x)^2$ ; the demo and part (b) above use $C=400$ . To the nearest 10, find a critical value $C_0$ such that the solution approaches zero asymptotically if $C < C_0$ , and infinity otherwise. (Use the solution at $t=5$ to determine the critical value.)

Exercise 11.5.4

⌨ The Allen–Cahn equation is used as a model for systems that prefer to be in one of two stable states. The governing PDE is

u_t = u(1-u^2) + \epsilon u_{xx}.

(11.5.13)

For this problem, assume $\epsilon=10^{-3}$ , $-1\le x \le 1$ , boundary conditions $u(\pm 1,t) = -1$ and initial condition

u(x,0) = -1 + \beta (1-x^2) e^{-20x^2},

(11.5.14)

where β is a parameter. Use Function 11.5.2 with $m=199$ .

(a) Solve the problem with $\beta=1.1$ up to time $t=8$ , plotting the solution at 6 equally spaced times. (You should see the solution decay down to the constant value -1.)

(b) Solve again with $\beta = 1.6$ . (This time the part of the bump will grow to just about reach $u=1$ , and stay there.)

Exercise 11.5.5

⌨ The Fisher equation is $u_t=u_{xx} + u - u^2$ . Assume that $0\le x \le 6$ and that the boundary conditions are $u_x(0,t)=u(6,t)=0$ .

(a) For the initial condition $u(x,0) = \frac{1}{2}[1+\cos(\pi x/2)]$ , use Function 11.5.2 with $m=80$ to solve the Fisher equation and plot the solution at times $t=0,0.5,\ldots,3$ . What is $u(0,3)$ ?

(b) Repeat part (a), but increase the final time until it appears that the solution has reached a steady state (i.e., stopped changing in time). Find an accurate value of $u(0,t)$ as $t \to \infty$ .

(c) If we set $u_t=0$ in the Fisher equation at steady state, we get a TPBVP in $u(x)$ . Use Function 10.5.1 with $m=300$ to solve this BVP, and make sure that the value at $x=0$ matches the result of part (b) to at least four digits.

Footnotes¶

An important advanced feature of Julia is multiple dispatch, which allows you to make multiple definitions of a function for different sequences and types of input arguments. Thus, addition to the original Function 11.5.2, we could also define a modified version in which g₁ and g₂ are of numeric type for the Dirichlet case. The correct version would be chosen (dispatched) depending on how the boundary conditions were supplied by the caller, allowing us speed when possible and generality as a fallback.
↩