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
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
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:
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:
The boundary conditions are used only in the definition of \(\mathbf{f}\). As in Section 10.5, define
Then Equation (11.5.2) takes the form
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.
Recall the Black–Scholes PDE (11.1.7),
subject to \(u(0)=0\) and \(u_x(S_\text{max})=1\). These imply \(u_0=0\) and
Hence
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
Hence
Implementation#
The steps to evaluate \(\mathbf{f}\) in (11.5.4) now go as follows.
Given a value of \(t\) and \(\mathbf{v}\),
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)\).
Solution of parabolic PDEs by the method of lines
1"""
2 parabolic(ϕ,xspan,m,g₁,g₂,tspan,init)
3
4Solve a parabolic PDE by the method of lines. The PDE is
5∂u/∂t = `ϕ`(t,x,u,∂u/∂x,∂^2u/∂x^2), `xspan` gives the space
6domain, m gives the degree of a Chebyshev spectral discretization,
7`g₁` and `g₂` are functions of (u,∂u/∂x) at the domain ends that
8should be made zero, `tspan` is the time domain, and `init` is a
9function of x that gives the initial condition. Returns a vector
10`x` and a function of t that gives the semidiscrete solution at `x`.
11"""
12function parabolic(ϕ,xspan,m,g₁,g₂,tspan,init)
13 x,Dₓ,Dₓₓ = diffcheb(m,xspan)
14 int = 2:m # indexes of interior nodes
15
16 function extend(v)
17 function objective(ubc)
18 u₀,uₘ = ubc
19 uₓ = Dₓ*[u₀;v;uₘ]
20 return [g₁(u₀,uₓ[1]),g₂(uₘ,uₓ[end])]
21 end
22 ubc = levenberg(objective,[0,0])[end]
23 return [ubc[1];v;ubc[2]]
24 end
25
26 function ode!(f,v,p,t)
27 u = extend(v)
28 uₓ,uₓₓ = Dₓ*u,Dₓₓ*u
29 @. f = ϕ(t,x[int],u[int],uₓ[int],uₓₓ[int])
30 end
31
32 ivp = ODEProblem(ode!,init.(x[int]),float.(tspan))
33 u = solve(ivp)
34
35 return x,t->extend(u(t))
36end
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
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
Show code cell source
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 /Users/driscoll/Documents/fnc-julia/diffusion/boundaries-heat.mp4
We solve a heat equation with a nonlinear source term,
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);
Show code cell source
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 /Users/driscoll/Documents/fnc-julia/diffusion/boundaries-source.mp4
Finally, we return to the example of the Black–Scholes equation from Section 11.1.
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₀);
Show code cell source
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 /Users/driscoll/Documents/fnc-julia/diffusion/boundaries-bs.mp4
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₁
andg₂
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.