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).

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:

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 Nonlinearity and boundary conditions, define

Then Equation (11.5.2) takes the form

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.

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)$.

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.2 is some speed in such special cases.^[1]

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

11.5.3Exercises¶

1. ✍ 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.6) for $u_0$ and $u_m$.

2. ⌨ Use Function 11.5.2 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)$.

3. Consider Demo 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=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.

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.)

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.

6. ⌨ Modify Function 11.5.2 for the special case of Dirichlet boundary conditions, in which (11.5.2) becomes simply

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

   (11.5.15)

   Your function should accept numbers α and β as input arguments in place of $g_1$ and $g_2$. Test your function on the problem in Demo 11.5.3.

7. ⌨ Modify Function 11.5.2 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.$