Two-point BVP - Fundamentals of Numerical Computation

The initial-value problems of Chapter 6 are characterized by an ordinary differential equation plus a value of the solution’s state at one value of the independent variable (typically, time).

In a boundary-value problem, the state is not entirely given at any point. Instead, partial information is given at multiple values of the independent variable. We will focus on the most common type.

10.1.1Definition¶

Definition 10.1.1 (Two-point boundary-value problem: TPBVP)

\begin{split} u''(x) &= \phi(x,u,u'), \qquad a \le x \le b,\\ g_1(u(a),u'(a)) &= 0,\\ g_2(u(b),u'(b)) &= 0. \end{split}

(10.1.1)

The functions $g_1$ and $g_2$ are called boundary conditions.

This TPBVP is said to be linear if the dependence of ϕ, $g_1$ , and $g_2$ on the solution $u(x)$ is linear. Specifically, this means that

\phi(x, u, u') = p(x)u' + q(x)u + r(x)

(10.1.2)

for some coefficient functions $p$ , $q$ , and $r$ , and

g_i(u, u') = \alpha_i u + \beta_i u' - \gamma_i, \quad i=1,2,

(10.1.3)

for some constants $\alpha_i$ , $\beta_i$ , and $\gamma_i$ .

An IVP for the same ODE as in (10.1.1) would specify values for both $u(a)$ and $u'(a)$ . Although this may look like a minor difference, the characters of IVPs and BVPs are substantially different. In a typical IVP, in which the independent variable is often time, the initial value determines everything about the future course of the solution. In a TPBVP, however, complete information is not given at any one location in the domain, and there may be more than one way (or no way) to satisfy the boundary conditions.

Certain special cases of the boundary conditions have their own nomenclature.

Example 10.1.1 (Multiple solutions in a BVP)

An ideal pendulum of length $L$ satisfies the ODE $\theta'' + \frac{g}{L}\sin \theta = 0$ , where $\theta(t)$ is the angle of the pendulum’s rod from straight downward, and $g$ is gravitational acceleration.

If we pull the pendulum bob 1 radian and release it from rest, then we have an IVP with the initial condition $\theta(0)=1$ , $\theta'(0)=0$ . Everything about the future trajectory of the pendulum is completely determined by that condition. But if instead we want to know how far to pull up the pendulum bob initially so that it is in the downward position 2 seconds later, then we have the Neumann boundary condition $\theta'(0)=0$ and the Dirichlet boundary condition $\theta(2)=0$ .

Clearly, one solution is $\theta(0)=0$ , and the pendulum never moves! But there is also an initial deflection such that $\theta(2)=0$ momentarily as the pendulum finishes its first downward swing. There could be still more solutions in which the pendulum completes one or more half-periods before passing again through the downward position at just the right moment.

While time can be the independent variable in a TPBVP, as in Example 10.1.1, it is often space, which has no intrinsic direction of information flow.

Example 10.1.2 (MEMS actuator)

A micro-electromechanical system (MEMS) can have an actuator consisting of two flat, parallel, disk-shaped surfaces, one at $z=0$ and the other at $z=1$ . The surface at $z=0$ is a rigid metal plate. The surface at $z=1$ is an elastic membrane fixed only at its boundary. When the surfaces are given different electric potentials, the membrane deflects in response to the induced electric field, and the field is also a position of the deflection. Assuming circular symmetry, one may derive the ordinary differential equation Pelesko & Driscoll (2006)

\frac{d^2w}{d r^2} + \frac{1}{r}\, \frac{d w}{d r} = \frac{\lambda}{w^2},

(10.1.4)

where $w(r)$ is the vertical position of the membrane at radius $r$ , and λ is proportional to the applied electric potential. We will use $0\le r\le 1$ . There are also the supplemental physical conditions

w'(0)=0, \qquad w(1)=1,

(10.1.5)

derived from the circular symmetry and fixing the edge of the membrane, respectively. This is a nonlinear boundary-value problem. The solution should respect $0<w(r)\le 1$ in order to be physically meaningful.

10.1.2Numerical solution¶

We can solve the TPBVP (10.1.1) by recasting the problem as a first-order system, as introduced in Section 6.3.2.

Example 10.1.3 (Solving a TPBVP numerically)

As a system, the MEMS problem from Example 10.1.2 uses $y_1=w$ , $y_2=w'$ to obtain

\begin{split} y_1' &= y_2, \\ y_2' &= \frac{\lambda}{y_1^2} - \frac{y_2}{r}. \end{split}

(10.1.6)

This has the form $\mathbf{y}' = \mathbf{f}(r,\mathbf{y})$ . In the original variable, the boundary conditions are $w'(0)=0$ and $w(1)=1$ . We have to encode these in the new variable as $y_2(0)=0$ and $y_1(1)-1 = 0$ .

Julia

MATLAB

Python

Example 10.1.3

To solve this problem, we have to define functions for the ODE and boundary conditions. The first, which is similar to what we do for solving IVPs, returns the computed values of $y_1'$ and $y_2'$ given values of the dependent and independent variable. We will use an in-place form of this function here, where the desired output is actually the first input argument. (Notice that the return value is irrelevant with the in-place style.)

function ode!(f, y, λ, r)
    f[1] = y[2]
    f[2] = λ / y[1]^2 - y[2] / r
    return nothing
end;

To encode the boundary conditions $y_2(0)=0$ , $y_1(2)-1=0$ , we define a function for their residual values. In other words, these are quantities that the solver has to make zero in order to satisfy the boundary conditions. This uses in-place style as well.

function bc!(g, y, λ, r)    # output, solution vector, parameter, indep. var.
    g[1] = y[1][2]          # first node, second component = 0
    g[2] = y[end][1] - 1    # last node, first component = 1
    return nothing
end;

In the bc! function, the y argument is just like an IVP solution from Basics of IVPs. Thus, y(0) is the value of the solution at $x=0$ , and the second component of that value is what we wish to make zero. Similarly, y(1)[1] is the notation for $y_1(1)$ , which is supposed to equal 1.

The domain of the mathematical problem is $r\in [0,1]$ . However, there is a division by $r$ in the ODE, so we want to avoid $r=0$ by truncating the domain a bit.

domain = (1e-15, 1.0)

(1.0e-15, 1.0)

We need one last ingredient that is not part of the mathematical specification: an initial guess for the solution $\mathbf{y}(r)$ . As we will see, this plays the same role as initialization in Newton’s method for rootfinding. Here, we try a constant value for each component.

est = [1, 0]

2-element Vector{Int64}:
 1
 0

Now we set up and solve a BVProblem with the parameter value $\lambda=0.6$ .

using BoundaryValueDiffEq, OrdinaryDiffEq, Plots
bvp = BVProblem(ode!, bc!, est, domain, 0.6)
y = solve(bvp, Shooting(Tsit5()))
plot(y;
    label = [L"w" L"w'"],
    legend = :right,
    xlabel = L"r",
    ylabel = "solution",
    title = "Solution of MEMS problem for λ=0.6",
)

To visual accuracy, the boundary conditions have been enforced. We can check them numerically.

@show y(0)[2];    # y vector at r=0, 2nd component
@show y(1)[1];    # y vector at r=1, 1st component

(y(0))[2] = -9.668656384782223e-16


(y(1))[1] = 0.9999999999999035

If the solver failed to satisfy these conditions, we would have to try a different initial guess.

Example 10.1.3

To define the BVP, we need to define some functions for the ODE and boundary conditions. (For this simple problem, we will use anonymous functions, but for a more substantial one, it would be better to use named functions.) The first, which is similar to what we do for solving IVPs, returns the computed values of $y_1'$ and $y_2'$ given values of the dependent and independent variable.

lambda = 0.6;
bvpfcn = @(r, y) [ y(2); lambda / y(1)^2 - y(2) / r  ];    % column vector [y1'(r); y2'(r)]

bcfcn = @(ya, yb) [ ya(2); yb(1) - 1 ];    % y_2(a) = 0;  y_1(b) = 1

We will define the domain next. Since there is a division by $r$ in the ODE, we want to avoid $r=0$ by truncating the domain a tiny bit. We use equispaced nodes here.

nodes = linspace(eps, 1, 61);

We need to make one more defintion that isn’t part of the mathematical formulation. Rather, it provides an initial guess for the solution. As we will see, this plays the same role as initialization in Newton’s method for rootfinding. Here we choose both components to be constant.

y_init = @(r) [ 1; 0 ];
sol_init = bvpinit(nodes, y_init);    % create the initial guess

Now we can solve the BVP using the bvp4c function.

sol = bvp4c(bvpfcn, bcfcn, sol_init);

plot(sol.x, sol.y, '-o')
xlabel('r'), ylabel('y(r)')
title('Solution of the membrane problem')
legend("w(r)", "w'(r)", location="east");

It’s smart to check visually that the boundary conditions are satisfied. The left BC implies a horizontal tangent at $r=0$ , and the right BC implies that $w(1)=1$ . If the solver failed to satisfy these conditions, we would have to try a different initial guess.

Example 10.1.3

lamb = 0.6
def ode(r, y):
    return array([
        y[1],
        lamb / y[0]**2 - y[1] / r
    ])

def bc(ya, yb):    # given vectors y(a), y(b)
    return array([
        ya[1],
        yb[0] - 1
    ])

The domain of the mathematical problem is $r\in [0,1]$ . However, there is a division by $r$ in this ODE, so we want to avoid $r=0$ by truncating the domain a bit.

a, b = finfo(float).eps, 1

r = linspace(a, b, 50)
w = ones(r.size)         # w(r) = 1 everywhere
dw_dr = zeros(r.size)    # w'(r) = 0 everywhere
y_init = vstack([w, dw_dr])

Now we can solve the problem using solve_bvp from scipy.integrate.

from scipy.integrate import solve_bvp
sol = solve_bvp(ode, bc, r, y_init)
print(f"Solved at {sol.x.size} nodes.")
plot(sol.x, sol.y[0])
xlabel("$r$"),  ylabel("$w(r)$")
title("Solution of MEMS problem for $\\lambda=0.6$");

Solved at 998 nodes.

It’s important to visually check the boundary conditions. The left BC implies a horizontal tangent at $r=0$ , and the right BC implies that $w(1)=1$ . If the solver failed to satisfy these conditions, we would have to try a different initial guess.

Characterizing the conditioning of a TPBVP theoretically is difficult. There are some numerical tools going by the name sensitivity analysis, but the details are too lengthy for us to explore here.

10.1.3Exercises¶

Exercise 10.1.2

✍ For each BVP, verify that the given solution is valid, i.e., check that the differential equation and the boundary conditions are satisfied.

(a) $u'' - 2xu' + 8 u = 0, \quad u(0) = 1, \; u(1) = -\frac{5}{3}\qquad$ Hermite equation

Solution: $u(x) = \frac{4}{3}x^4-4x^2+1$

(b) $xu'' + (1-x) u' + 3 u = 0, \quad u(0) =1, \; u'(1) =-\frac{1}{2} \qquad$ Laguerre equation

Solution: $u(x) = \frac{1}{6}(6-18x+9x^2-x^3)$

(c) $(1-x^2)u'' - xu' + 25u = 0, \quad u'(0) = 5, \; u\left(\frac{1}{2}\right) = \frac{1}{2} \qquad$ Chebyshev equation

Solution: $u(x) = T_5(x) = 16x^5-20x^3+5x$

(d) $(1-x^2) u'' -2xu' + 12 u = 0, \quad u(0)+2u'(0) = -3, \; u(1) = 1 \qquad$ Legendre equation

Solution: $u(x) = P_3(x) = \frac{1}{2}(5x^3-3x)$

(e) $u u'' = 3(u')^2, \quad u(-1) = 1, \; u(2) = \frac{1}{2}$

Solution: $u(x) = ( x+2 )^{-1/2}$

Exercise 10.1.5

⌨ The stationary Allen–Cahn equation is a model of phase changes, such as the change from liquid to solid. In one spatial dimension it can be written as

\epsilon u'' = u^3-u, \qquad 0 \le x \le 1, \qquad u(0)=-1, \quad u(1)=1.

(10.1.7)

As $\epsilon\rightarrow 0$ , the solution tends toward a step function transition between -1 and 1. By symmetry, $u'(x)=-u'(1-x)$ .

(a) Use Function 10.2.1 with initial solution estimates $(u=-1,u'=0)$ to solve the equation for $\epsilon=0.2$ . Plot the solution.

(b) Repeat part (a) for $\epsilon=0.02$ .

(c) Repeat part (a) for $\epsilon=0.002$ . (This is a difficult problem for the default method.) Try a few different initializations for $u$ , and plot all the results. Do any seem to be valid solutions of the BVP?

References¶

Pelesko, J. A., & Driscoll, T. A. (2006). The Effect of the Small-Aspect-Ratio Approximation on Canonical Electrostatic MEMS Models. Journal of Engineering Mathematics, 53(3–4), 239–252. 10.1007/s10665-005-9013-2
Carrier, G. F. (1970). Singular Perturbation Theory and Geophysics. SIAM Review, 12(2), 175–193. 10.1137/1012041

Preface

10. Boundary-value problems

Preface

Shooting