from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

10.5.Nonlinearity and boundary conditions¶

Collocation for nonlinear differential equations operates on the same principle as for linear problems: replace functions by vectors and replace derivatives by differentiation matrices. But because the differential equation is nonlinear, the resulting algebraic equations are as well. We will therefore need to use a quasi-Newton or similar method as part of the solution process.

We consider the TPBVP (10.1.1), reproduced here:

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

As in Collocation for linear problems, the function $u(x)$ is replaced by a vector $\mathbf{u}$ of its approximated values at nodes $x_0,x_1,\ldots,x_n$ (see Equation (10.4.2)). We define derivatives of the sampled function as in (10.4.3) and (10.4.4), using suitable differentiation matrices $\mathbf{D}_x$ and $\mathbf{D}_{xx}$ .

The collocation equations, ignoring boundary conditions for now, are

\mathbf{D}_{xx} \mathbf{u} - \mathbf{r}(\mathbf{u}) = \boldsymbol{0},

(10.5.2)

where

r_i(\mathbf{u}) = \phi(x_i,u_i,u_i'), \qquad i=0,\ldots,n.

(10.5.3)

and $\mathbf{u}'=\mathbf{D}_x\mathbf{u}$ .

We impose the boundary conditions in much the same way as in Collocation for linear problems. Again define the rectangular boundary removal matrix $\mathbf{E}$ as in (10.4.8), and replace the equations in those two rows by the boundary conditions:

\mathbf{f}(\mathbf{u}) = \begin{bmatrix} \mathbf{E} \bigl( \mathbf{D}_{xx}\mathbf{u} - \mathbf{r}(\mathbf{u}) \bigr) \\[1mm] g_1(u_0,u_0') \\[1mm] g_2(u_n,u_n') \end{bmatrix} = \boldsymbol{0}.

(10.5.4)

The left-hand side of (10.5.4) is a nonlinear function of the unknowns in the vector $\mathbf{u}$ , so (10.5.4) is an $(n+1)\times 1$ set of nonlinear equations, amenable to solution by the techniques of Chapter 4.

Example 10.5.1

Given the BVP

u'' - \sin(xu) + \exp(xu')=0, \quad u(0)=-2, \; u'(3/2)=1,

(10.5.5)

we compare to the standard form (10.5.1) and recognize

\phi(x,u,u') = \sin(xu)-\exp(xu').

(10.5.6)

Suppose $n=3$ for an equispaced grid, so that $h=\frac{1}{2}$ , $x_0=0$ , $x_1=\frac{1}{2}$ , $x_2=1$ , and $x_3=\frac{3}{2}$ . There are four unknowns. We compute

\begin{gather*} \mathbf{D}_{xx} = \frac{1}{1/4} \begin{bmatrix} 2 & -5 & 4 & -1 \\ 1 & -2 & 1 & 0 \\ 0 & 1 & -2 & 1 \\ -1 & 4 & -5 & 2 \end{bmatrix}, \quad \mathbf{D}_x = \frac{1}{1} \begin{bmatrix} -3 & 4 & -1 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & -1 & 0 & 1 \\ 0 & 1 & -4 & 3 \end{bmatrix}, \\[3mm] \mathbf{E} \mathbf{r}(\mathbf{u}) = \begin{bmatrix} \sin\left(\frac{u_1}{2}\right) - \exp\left(\frac{u_2-u_0}{2}\right) \\[1mm] \sin(u_2) - \exp\left( u_3-u_1 \right) \end{bmatrix}, \\[3mm] \mathbf{f}(\mathbf{u}) = \begin{bmatrix} (4u_0 -8u_1 + 4u_2) - \sin\left(\frac{u_1}{2}\right) + \exp\left(\frac{u_2-u_0}{2}\right) \\[1mm] (4u_1 -8u_2 + 4u_3) - \sin(u_2) + \exp\left( u_3-u_1 \right) \\[1mm] u_0 + 2 \\[1mm] (u_1 - 4u_2 + 3u_3) - 1 \end{bmatrix}. \end{gather*}

10.5.1Implementation¶

Our implementation using second-order finite differences is Function 10.5.1. It’s surprisingly short, considering how general it is, because we have laid a lot of groundwork already.

Algorithm 10.5.1 (bvp)

bvp.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
def bvp(phi, xspan, ga, gb, init):
    """
    bvp(phi, xspan, ga, gb, init)

    Use finite differences to solve a two-point boundary value problem. 
    The ODE is u'' = phi(x, u, u') for x in (a,b). The functions 
    ga(u(a), u'(a)) and gb(u(b), u'(b)) specify the boundary conditions. 
    The value init is an initial guess for [u(a), u'(a)].

    Return vectors for the nodes and the values of u.
    """
    n = len(init) - 1
    x, Dx, Dxx = diffmat2(n, xspan)
    h = x[1] - x[0]
    def residual(u):
        # Compute the difference between u'' and phi(x,u,u') at the
        # interior nodes and appends the error at the boundaries.
        du_dx = Dx @ u  # discrete u'
        d2u_dx2 = Dxx @ u  # discrete u''
        f = d2u_dx2 - phi(x, u, du_dx)

        # Replace first and last values by boundary conditions.
        f[0] = ga(u[0], du_dx[0]) / h
        f[n] = gb(u[n], du_dx[n]) / h
        return f

    u = levenberg(residual, init.copy())
    return x, u[-1]

About the code

The nested function residual uses differentiation matrices computed externally to it, rather than computing them anew on each invocation. As in Function 10.4.1, there is no need to form the row-deletion matrix $\mathbf{E}$ explicitly. In lines 23--24, we divide the values of $g_1$ and $g_2$ by a factor of $h$ . This helps scale the residual components more uniformly and improves the robustness of convergence a bit.

In order to solve a particular problem, we must write a function that computes $\phi$ for vector-valued inputs $\mathbf{x}$ , $\mathbf{u}$ , and $\mathbf{u}'$ , and functions for the boundary conditions. We also have to supply init, which is an estimate of the solution used to initialize the quasi-Newton iteration. Since this argument is a vector of length $n+1$ , it sets the value of $n$ in the discretization.

Example 10.5.2 (BVP for a nonlinear pendulum)

Suppose a damped pendulum satisfies the nonlinear equation $\theta'' + 0.05\theta'+\sin \theta =0$ . We want to start the pendulum at $\theta=2.5$ and give it the right initial velocity so that it reaches $\theta=-2$ at exactly $t=5$ . This is a boundary-value problem with Dirichlet conditions $\theta(0)=2.5$ and $\theta(5)=-2$ .

The first step is to define the function $\phi$ that equals $\theta''$ .

phi = lambda t, theta, omega: -0.05 * omega - sin(theta)

Next, we define the boundary conditions.

ga = lambda u, du: u - 2.5
gb = lambda u, du: u + 2

The last ingredient is an initial estimate of the solution. Here we choose $n=100$ and a linear function between the endpoint values.

init = linspace(2.5, -2, 101)

We find a solution with negative initial slope, i.e., the pendulum is initially pushed back toward equilibrium.

t, theta = FNC.bvp(phi, [0, 5], ga, gb, init)
plot(t, theta)
xlabel("$t$")
ylabel("$\theta(t)$")
title("Pendulum over [0,5]");

If we extend the time interval longer for the same boundary values, then the initial slope must adjust.

t, theta = FNC.bvp(phi, [0, 8], ga, gb, init)
plot(t, theta)
xlabel("$t$")
ylabel("$\theta(t)$")
title("Pendulum over [0,8]");

This time, the pendulum is initially pushed toward the unstable equilibrium in the upright vertical position before gravity pulls it back down.

The initial solution estimate can strongly influence how quickly a solution is found, or whether the quasi-Newton iteration converges at all. In situations where multiple solutions exist, the initialization can determine which is found.

Example 10.5.3 (BVP for a nonlinear MEMS device)

We look for a solution to the parameterized membrane deflection problem from Example 10.1.2,

w''+ \frac{1}{r}w'= \frac{\lambda}{w^2},\quad w'(0)=0,\; w(1)=1.

(10.5.7)

Here is the problem definition. We use a truncated domain to avoid division by zero at $r=0$ .

lamb = 0.5
phi = lambda r, w, dwdr: lamb / w**2 - dwdr / r
a, b = finfo(float).eps, 1
ga = lambda w, dw: dw
gb = lambda w, dw: w - 1

First we try a constant function as the initialization.

init = ones(201)
r, w1 = FNC.bvp(phi, [a, b], ga, gb, init)
plot(r, w1)
fig, ax = gcf(), gca()
xlabel("$r$"),  ylabel("$w(r)$")
title("Solution of the MEMS problem");

It’s not necessary that the initialization satisfy the boundary conditions. In fact, by choosing a different constant function as the initial guess, we arrive at another valid solution.

r, w2 = FNC.bvp(phi, [a, b], ga, gb, 0.5 * init)
ax.plot(r, w2)
ax.set_title("Multiple solutions of the MEMS problem");
fig

10.5.2Parameter continuation¶

Sometimes the best way to get a useful initialization is to use the solution of a related easier problem, a technique known as parameter continuation. In this approach, one solves the problem at an easy parameter value, and gradually changes the parameter value to the desired value. After each change, the most recent solution is used to initialize the iteration at the new parameter value.

Example 10.5.4 (Allen–Cahn equation)

We solve the stationary Allen–Cahn equation,

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

(10.5.8)

phi = lambda x, u, dudx: (u**3 - u) / epsilon
ga = lambda u, du: du
gb = lambda u, du: u - 1

Finding a solution is easy at larger values of $\epsilon$ .

epsilon = 0.05
init = linspace(-1, 1, 141)
x, u1 = FNC.bvp(phi, [0, 1], ga, gb, init)

plot(x, u1, label="$\\epsilon = 0.05$")
fig, ax = gcf(), gca()
xlabel("$x$"),  ylabel("$u(x)$")
legend(),  title("Allen-Cahn solution");

/Users/driscoll/miniforge3/envs/myst/lib/python3.13/site-packages/fncbook/chapter04.py:167: UserWarning: Iteration did not find a root.
  warnings.warn("Iteration did not find a root.")

Finding a good initialization is not trivial for smaller values of $\epsilon$ . But the iteration succeeds if we use the first solution as the initialization at the smaller $\epsilon$ .

epsilon = 0.002
x, u2 = FNC.bvp(phi, [0, 1], ga, gb, u1)
ax.plot(x, u2, label="$\\epsilon = 0.002$")
ax.legend()
fig

In this case we can continue further.

ϵ = 0.0005
x, u3 = FNC.bvp(phi, [0, 1], ga, gb, u2)
ax.plot(x, u3, label="$\\epsilon = 0.005$")
ax.legend()
fig

10.5.3Exercises¶

Exercise 10.5.7

⌨ The following nonlinear BVP was proposed by Carrier (for the special case $b=1$ in Carrier (1970)):

\epsilon u'' + 2(1-x^2)u +u^2 = 1, \quad u(-1) = u(1) = 0.

(10.5.12)

In order to balance the different components of the residual, it’s best to implement each boundary condition numerically as $u/\epsilon=0$ .

(a) Use Function 10.5.1 to solve the problem with $\epsilon=0.003$ , $n=200$ , and an initial estimate of all zeros. Plot the result; you should get a solution with 9 local maxima.

(b) Starting with the result from part (a) as an initialization, continue the parameter through the sequence

\epsilon = 3\times 10^{-3}, 3\times 10^{-2.8}, 3\times 10^{-2.6},\ldots, 3\times 10^{-1}.

(10.5.13)

The most recent solution should be used as the initialization for each new value of $\epsilon$ . Plot the end result for $\epsilon=0.3$ ; it should have one interior local maximum.

(c) Starting with the last solution of part (b), reverse the continuation steps to return to $\epsilon=0.003$ . Plot the result, which is an entirely different solution from part (a).

References¶

Ascher, U. M., & Petzold, L. R. (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Society for Industrial and Applied Mathematics. 10.1137/1.9781611971392
Carrier, G. F. (1970). Singular Perturbation Theory and Geophysics. SIAM Review, 12(2), 175–193. 10.1137/1012041