11.5. Boundaries#

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

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

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

(11.5.2)#\[\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}\]

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

Boundary 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:

(11.5.3)#\[\begin{split} \mathbf{u} = \begin{bmatrix} u_0 \\ \mathbf{v} \\ u_m \end{bmatrix},\end{split}\]

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:

(11.5.4)#\[ \frac{d \mathbf{v}}{d t} = \mathbf{f}(t,\mathbf{v}).\]

The boundary conditions are used only in the definition of \(\mathbf{f}\). As in Section 10.5, define

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

Then Equation (11.5.2) takes the form

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

Given a value of \(\mathbf{v}\) for the interior nodes, Equation (11.5.5) 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, \]

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. \]

Hence

\[\begin{split} \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}. \end{split}\]

Example 11.5.2

Returning to Example 11.5.1, suppose we use a global Chebyshev differentiation matrix for \(\mathbf{D}_x\) in (11.5.5). Then \(u_0=0\) and

\[ D_{m0} u_0 + D_{m1}u_1 + \cdots + D_{mm}u_m = 1. \]

Hence

\[\begin{split} \begin{bmatrix} 1 & 0 \\ D_{m0} & D_{mm} \end{bmatrix} \begin{bmatrix} u_0 \\ u_m \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} - \begin{bmatrix} 0 & \cdots & 0 & 0 \\ D_{m1} & \cdots & D_{m,m-2}& D_{m,m-1} \end{bmatrix} \mathbf{v}. \end{split}\]

Implementation#

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

Algorithm 11.5.3 : Time derivative for parabolic PDE

Given a value of \(t\) and \(\mathbf{v}\),

Use (11.5.5) to solve for \(u_0\) and \(u_m\).
Assemble the total vector \(\mathbf{u}\) from (11.5.3).
Use the spatial semidiscretization to evaluate \(\phi\) at all the nodes.
Ignore the boundary nodes to get the value of \(\mathbf{f}(t,\mathbf{v})\) in (11.5.4).

Our full implementation of the method of lines for (11.5.1)–(11.5.2) is given in Function 11.5.4. It uses Function 10.3.3 (diffcheb) to set up a Chebyshev discretization. The nested function extend performs steps 1–2 of Algorithm 11.5.3 by calling Function 4.6.4 (levenberg) to solve the potentially nonlinear system (11.5.5). Then it sets up and solves an IVP, adding steps 3–4 of Algorithm 11.5.3 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)\).

Function 11.5.4 : parabolic

Solution of parabolic PDEs by the method of lines

"""
    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

About the code

Line 29 uses the macro @. to assign into the vector f elementwise. Without it, the function would allocate space for the result of phi and then change f to point at that vector, and that would defeat the purpose of using the preallocated f for speed.

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 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.4 is some speed in such special cases.1

Demo 11.5.5

Let’s solve the heat equation on \([-1,1]\). First, we define the PDE \(u_t=u_{xx}\) and Dirichlet boundary conditions \(u(-1,t)=0\), \(u(1,t)=2\).

ϕ = (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 + sin(π*x/2) + 3*(1-x^2)*exp(-4x^2);

Now we can use Function 11.5.4 to solve the problem.

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

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

┌ Info: Saved animation to 
│   fn = /Users/driscoll/repos/fnc-julia/diffusion/boundaries-heat.mp4
└ @ Plots /Users/driscoll/.julia/packages/Plots/8K4be/src/animation.jl:114

Demo 11.5.6

We solve a heat equation with a nonlinear source term,

\[ u_t = u^2 + u_{xx}. \]

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.

ϕ = (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);

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

┌ Info: Saved animation to 
│   fn = /Users/driscoll/repos/fnc-julia/diffusion/boundaries-source.mp4
└ @ Plots /Users/driscoll/.julia/packages/Plots/8K4be/src/animation.jl:114

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

Demo 11.5.7

Here is the Black–Scholes PDE, with a homogeneous Dirichlet condition at zero stock price and a nonhomogeneous Neumann condition at the truncated right boundary.

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₀);

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

┌ Info: Saved animation to 
│   fn = /Users/driscoll/repos/fnc-julia/diffusion/boundaries-bs.mp4
└ @ Plots /Users/driscoll/.julia/packages/Plots/8K4be/src/animation.jl:114

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.

Exercises#

✍ Suppose second-order finite differences with \(m=3\) are used to discretize the heat equation on \(x \in [0,2]\), with boundary conditions \(u_x(0,t)=0\) and \(u_x(2,t)=1\). If at some time \(t\), \(u_1=1\) and \(u_2=2\), set up and solve the equations (11.5.5) for \(u_0\) and \(u_m\).
⌨ Use Function 11.5.4 to solve the heat equation for \(0\le x \le 5\) with initial condition \(u(x,0)=x(5-x)\) and subject to the boundary conditions \(u(0,t)=0\), \(u(5,t)-u_x(5,t)=5\). Plot the solution at \(t=1\) and find the value of \(u(2.5,1)\).
Consider Demo 11.5.6, 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=1\). What happens?

(c) ⌨ Let the initial condition be \(u(x,0) = C x^4(1-x)^2\); the demo uses \(C=400\). To the nearest 10, find a critical value \(C_0\) such that the solution approaches zero asymptotically if \(C < C_0\), but not otherwise.
⌨ 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}.\]

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},\]

where \(\beta\) is a parameter. Use Function 11.5.4 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.)
⌨ 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.4 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.2 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.
⌨ Modify Function 11.5.4 for the special case of Dirichlet boundary conditions, in which (11.5.2) becomes simply

\[ u(a,t) = \alpha,\; u(b,t)=\beta. \]

Your function should accept numbers \(\alpha\) and \(\beta\) as input arguments in place of \(g_1\) and \(g_2\). Test your function on the problem in Demo 11.5.5.
⌨ Modify Function 11.5.4 for the special case of homogeneous Neumann boundary conditions. There is no longer any need for the input arguments for \(g_1\) and \(g_2\). Your implementation should solve a \(2\times 2\) linear system of equations to find the boundary values within the nested function extend. Test your function on the heat equation on \(x \in [0,4]\), \(t\in [0,1]\) with initial condition \(u(x,0)=x^2(4-x)^4.\)

1: 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.4, 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.

Fundamentals of Numerical Computation

Boundaries

Contents

11.5. Boundaries#

Boundary removal#

Implementation#

Exercises#