Nonlinearity and boundary conditions - Fundamentals of Numerical Computation

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}, \\[2mm] \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}, \\[2mm] \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.7)

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)

Julia

MATLAB

Python

Solution of a nonlinear boundary-value problem

bvp.jl

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
"""
    bvplin(p, q, r, xspan, lval, rval, n)

Use finite differences to solve a linear bopundary value problem.
The ODE is u''+`p`(x)u'+`q`(x)u = `r`(x) on the interval `xspan`,
with endpoint function values given as `lval` and `rval`. There will
be `n`+1 equally spaced nodes, including the endpoints.

Returns vectors of the nodes and the solution values.
"""
function bvplin(p, q, r, xspan, lval, rval, n)
    x, Dₓ, Dₓₓ = diffmat2(n, xspan)

    P = diagm(p.(x))
    Q = diagm(q.(x))
    L = Dₓₓ + P * Dₓ + Q     # ODE expressed at the nodes

    # Replace first and last rows using boundary conditions.
    z = zeros(1, n)
    A = [[1 z]; L[2:n, :]; [z 1]]
    b = [lval; r.(x[2:n]); rval]

    # Solve the system.
    u = A \ b
    return x, u
end

Solution of a nonlinear boundary-value problem

bvp.m

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
29
30
31
32
33
34
  function [x,u] = bvp(phi, a, b, ga, gb, init)
%BVP      Solve a boundary-value problem by finite differences
%         with either Dirichlet or Neumann BCs.
% Input:
%   phi      defines u'' = phi(x,u,u') (function)
%   a, b     endpoints of the domain (scalars)
%   ga       residual boundary function of u(a), u'(a) 
%   gb       residual boundary function of u(b), u'(b) 
%   init     initial guess for the solution (length n+1 vector)
% Output:
%   x        nodes in x (vector, length n+1)
%   u        values of u(x)  (vector, length n+1)
%   res      function for computing the residual

n = length(init) - 1;
[x, Dx, Dxx] = diffmat2(n, [a, b]);
h = x(2) - x(1);

u = levenberg(@residual, init);
u = u(:, end);

    function f = residual(u)
        % Computes 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(1) =   ga(  u(1),   du_dx(1)) / h;
        f(end) = gb(u(end), du_dx(end)) / h;
    end

end

Solution of a nonlinear boundary-value problem

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]

In order to solve a particular problem, we must write a function that computes ϕ 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$ .

Julia

MATLAB

Python

Example 10.5.2

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

ϕ = (t, θ, ω) -> -0.05 * ω - sin(θ);

Next, we define the boundary conditions.

g₁(u, du) = u - 2.5
g₂(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.

The collect function turns a range object into a true vector.

init = collect(range(2.5, -2, length = 101));

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

using Plots
t, θ = FNC.bvp(ϕ, [0, 5], g₁, g₂, init)
plot(t, θ;
    xaxis=(L"t"),  yaxis=(L"\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, θ = FNC.bvp(ϕ, [0, 8], g₁, g₂, init)
plot(t, θ;
    xaxis=(L"t"),  yaxis=(L"\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.

Example 10.5.2

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

phi = @(t,theta,omega) -0.05 * omega - sin(theta);

Next, we define the boundary conditions.

ga = @(u, du) u - 2.5;
gb = @(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] = bvp(phi, 0, 5, ga, gb, init);
clf,  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] = 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.

Example 10.5.2

The first step is to define the function ϕ 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.8)

Julia

MATLAB

Python

Example 10.5.3

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

domain = [eps(), 1]
λ = 0.5
ϕ = (r, w, dwdr) -> λ / w^2 - dwdr / r
g₁(w, dw) = dw
g₂(w, dw) = w - 1;

First we try a constant function as the initialization.

init = ones(301)
r, w₁ = FNC.bvp(ϕ, domain, g₁, g₂, init)

plot(r, w₁;
    xaxis = (L"r"),  yaxis = (L"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.

init = 0.5 * ones(301)
r, w₂ = FNC.bvp(ϕ, domain, g₁, g₂, init)
plot!(r, w₂, title = "Two solutions of the MEMS problem")

Example 10.5.3

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

lambda = 0.5;
phi = @(r,w,dwdr) lambda./w.^2 - dwdr./r;
ga = @(w, dw) dw;
gb = @(w, dw) w - 1;
a = eps;  b = 1;

First we try a constant function as the initialization.

init = ones(301, 1);
[r, w1] = bvp(phi, a, b, ga, gb, init);

clf,  plot(r, w1)
xlabel('r'),  ylabel('w(r)')
title('Solution of the MEMS BVP')

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.

init = 0.5 * ones(301, 1);
[r, w2] = bvp(phi, a, b, ga, gb, init);
hold on,  plot(r, w2)
title("Two solutions of the MEMS BVP")

Example 10.5.3

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

Julia

MATLAB

Python

Example 10.5.4

ϕ = (x, u, dudx) -> (u^3 - u) / ϵ;
g₁(u, du) = du
g₂(u, du) = u - 1;

Finding a solution is easy at larger values of ε.

ϵ = 0.05
init = collect(range(-1, 1, length = 141))
x, u₁ = FNC.bvp(ϕ, [0, 1], g₁, g₂, init)

plot(x, u₁;
    label=L"\epsilon = 0.05",  legend=:bottomright,
    xaxis=(L"x"),  yaxis=(L"u(x)"),
    title = "Allen–Cahn solution")

However, finding a good initialization is not trivial for smaller values of ε. Note below that the iteration stops without converging to a solution.

ϵ = 0.002;
x, z = FNC.bvp(ϕ, [0, 1], g₁, g₂, init);

┌ Warning: Maximum number of iterations reached.
└ @ FNCFunctions ~/Documents/FNCFunctions/src/chapter04.jl:176

The iteration succeeds if we use the first solution instead as the initialization here.

x, u₂ = FNC.bvp(ϕ, [0, 1], g₁, g₂, u₁)
plot!(x, u₂; label = L"\epsilon = 0.002")

In this case we can continue further.

ϵ = 0.0005
x, u₃ = FNC.bvp(ϕ, [0, 1], g₁, g₂, u₂)
plot!(x, u₃, label = L"\epsilon = 0.0005")

Example 10.5.4

epsilon = 0.05;
phi = @(x, u, du_dx) (u.^3 - u) / epsilon;
ga = @(u, du) du;
gb = @(u, du) u - 1;

Finding a solution is easy at larger values of ε.

init = linspace(-1, 1, 141)';
[x, u1] = bvp(phi, 0, 1, ga, gb, init);
clf,  plot(x, u1, displayname="\epsilon = 0.05")
xlabel('x'),  ylabel('u(x)')
title('Allen-Cahn solution') 
legend(location="northwest")

However, finding a good initialization is not trivial for smaller values of ε. Note below that the iteration stops without converging to a solution.

epsilon = 0.002;
phi = @(x, u, du_dx) (u.^3 - u) / epsilon;
[x, z] = bvp(phi, 0, 1, ga, gb, init);

The iteration succeeds if we use the first solution instead as the initialization here.

[x, u2] = bvp(phi, 0, 1, ga, gb, u1);
hold on,  plot(x, u2, displayname="\epsilon = 0.002")

In this case we can continue further.

epsilon = 0.0005;
phi = @(x, u, du_dx) (u.^3 - u) / epsilon;
[x, u3] = bvp(phi, 0, 1, ga, gb, u2);
plot(x, u3, displayname="\epsilon = 0.0005")

Example 10.5.4

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 = 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");

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

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¶

✍ This exercise is about the nonlinear boundary-value problem
$u'' = \frac{3(u')^2}{u} , \quad u(-1) = 1, \; u(2) = \frac{1}{2}.$
(10.5.10)
(a) Verify that the exact solution is $u(x) = ( x+2 )^{-1/2}$ .
(b) Write out the finite-difference approximation (10.5.4) with a single interior point ( $n=2$ ).
(c) Solve the equation of part (b) for the lone interior value $u_1$ .
⌨ (a) Use Function 10.5.1 to solve the problem of Exercise 1 for $n=80$ . In a 2-by-1 subplot array, plot the finite-difference solution and its error.
(b) ⌨ For each $n=10,20,40,\ldots,640$ , find the infinity norm of the error on the same problem. Make a log-log plot of error versus $n$ and include a graphical comparison to second-order convergence.
⌨ (Adapted from Ascher & Petzold (1998).) Use Function 10.5.1 twice with $n=200$ to solve
$u'' + e^{u+0.5} = 0, \quad y(0) = y(1) = 0,$
(10.5.11)
with initializations $7 \sin(x)$ and $\frac{1}{4} \sin(x)$ . Plot the solutions together on one graph.
⌨ Use Function 10.5.1 to compute the solution to the Allen–Cahn equation in Demo 10.5.4 with $\epsilon=0.02$ . Determine numerically whether it is antisymmetric around the line $x=0.5$ ---that is, whether $u(1-x)=-u(x)$ . You should supply evidence that your answer is independent of $n$ .
⌨ Consider the pendulum problem from Example 10.1.1 with $g=L=1$ . Suppose we want to release the pendulum from rest such that $\theta(5)=\pi/2$ . Use Function 10.5.1 with $n=200$ to find one solution that passes through $\theta=0$ , and another solution that does not. Plot $\theta(t)$ for both cases together.
⌨ The BVP
$u'' = x \operatorname{sign}(1-x) u, \quad u(-6)=1, \; u'(6)=0,$
(10.5.12)
forces $u''$ to be discontinuous at $x=1$ , so finite differences may not converge to the solution at their nominal order of accuracy.
(a) Solve the problem using Function 10.5.1 with $n=1400$ , and make a plot of the solution. Store the value at $x=6$ for use as a reference high-accuracy solution.
(b) For each $n=100,200,300,\ldots,1000$ , apply Function 10.5.1, and compute the error at $x=6$ . Compare the convergence graphically to second order.
⌨ 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.13)
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.14)
The most recent solution should be used as the initialization for each new value of ε. 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).
⌨ Demo 10.5.3 finds two solutions at $\lambda=0.5$ . Continue both solutions by taking 50 steps from $\lambda=0.5$ to $\lambda=0.79$ . Make a plot with λ on the horizontal axis and $w(0)$ on the vertical axis, with one point to represent each solution found. You should get two paths that converge as λ approaches 0.79 from below.

References¶

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

Preface

Collocation for linear problems

Preface

The Galerkin method