Two-point BVP - Fundamentals of Numerical Computation

The initial-value problems of 6. Initial-value problems for ODEs are characterized by an ordinary differential equation plus a value of the solution’s state at one value of the independent variable.

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.

Often the domain of $x$ in (10.1.1) is not explicitly stated but is implied by the definitions of $g_1$ and $g_2$ These functions are called boundary conditions. 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 change, IVPs and BVPs are quite different. Conceptually, the difference is like that between time and space. In a typical IVP, in which the independent variable is often time, the initial value determines everything about the future course of the solution. Even in a simple BVP, however, the necessary information is spread across the domain, and there may be more than one way to satisfy the boundary conditions.

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

Example 10.1.1

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

A micromechanical electrically driven actuator consists of two flat disk-shaped surfaces in parallel, 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.3)

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

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.1Numerical solution¶

We can solve the TPBVP (10.1.1) by recasting the problem as a first-order system in the usual way.

Example 10.1.3 (Solving a TPBVP)

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

This has the form $\mathbf{y}' = \mathbf{f}(r,\mathbf{y})$ . The boundary conditions are $y_2(0)=0$ and $y_1(1)=1$ .

Julia

MATLAB

Python

Example 10.1.3

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}

(1)

We will code an in-place form of this ODE, in which the first argument is used to return the computed values of $y_1'$ and $y_2'$ .

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

Notice that the return value is irrelevant with the in-place style. We use the same style for the boundary conditions $y_2(0)=0$ , $y_1(1)=1$ .

function bc!(g, y, λ, r)
    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 = (eps(), 1.0)

(2.220446049250313e-16, 1.0)

We need one last ingredient that is not part of the mathematical setup: an initial estimate for the solution. 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_2(0)
@show y(1)[1];    # y_1(1)

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


(y(1))[1] = 0.9999999999998375

Example 10.1.3

To define the BVP, we need to define some functions. (For this simple problem, we will use anonymous functions, but for a more substantial one, it would be better to use separate files.) The first defines $\mathbf{f}(x, \mathbf{y})$ .

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

Our second function defines the boundary conditions. It takes $\mathbf{y}(a)$ and $\mathbf{y}(b)$ as arguments and returns a vector of residuals; i.e., values that should be zero when the boundary conditions are satisfied.

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

The third function we define isn’t part of the mathematical formulation. Rather, it provides an initial guess for the solution. Here we choose both components to be constant.

y_init = @(r) [ 1; 0 ];

Now we can solve the BVP using the bvp4c function. We need to specify the nodes on which to solve the problem. The domain of the mathematical problem is $r\in [0,1]$ .But since there is a division by $r$ in the ODE, we want to avoid $r=0$ by truncating the domain a tiny bit.

nodes = linspace(eps, 1, 61);
sol_init = bvpinit(nodes, y_init);
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.

Example 10.1.3

To solve this problem, we have to define functions for the ODE and boundary conditions. The first returns the computed values of $y_1'$ and $y_2'$ .

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

To encode the boundary conditions $y_2(0)=0$ , $y_1(2)=1$ , we define a function for their residual values.

def bc(ya, yb):    # given 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 the 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)
y_init = vstack([ones(r.size), zeros(r.size)])

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.

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.2Exercises¶

✍ In each case, explain whether the TPBVP is linear or nonlinear.
(a) $x^2 u'' +xu' + (x^2 - 1) u = 0, \quad u(0) =1,\; u(4) =0 \qquad$ (Bessel equation)
(b) $u'' - u u' = 1, \quad u(0) = u'(1) = 1$
(c) $u'' + u = 1, \quad u(0) = u(1) u'(1) = 1$
(d) $\epsilon u'' +2(1-x^2) u +u^2 = 1, \quad u(-1) = u(1) = 0 \qquad$ (Carrier equation Carrier (1970))
(e) $u u'' = 3(u')^2, \quad u(-1) = 1, \; u(2) = \frac{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}$
⌨ For each TPBVP in Exercise 2, use solve to find the solution. Plot the solution and separately plot the error as a function of $x$ . (In some cases you will need to truncate the domain to avoid division by zero.)

⌨ Consider the pendulum from Example 10.1.1 with $g=L=1$ . Suppose we want to release the pendulum from rest such that $\theta(5)=\pi/2$ . Find one solution that passes through $\theta=0$ and another solution that does not. Plot $\theta(t)$ for both cases together. (Hint: Vary the initial estimate for the solution.)

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