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}, \\[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)

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
27
28
29
30
31
"""
    bvp(ϕ, xspan, g₁, g₂, init)

Finite differences to solve a two-point boundary value problem with
ODE u'' = `ϕ`(x,u,u') for x in `xspan`, left boundary condition
`g₁`(u,u')=0, and right boundary condition `g₂`(u,u')=0. The value
`init` is an initial estimate for the values of the solution u at
equally spaced values of x, which also sets the number of nodes.

Returns vectors for the nodes and the values of u.
"""
function bvp(ϕ, xspan, g₁, g₂, init)
    n = length(init) - 1
    x, Dₓ, Dₓₓ = diffmat2(n, xspan)
    h = x[2] - x[1]

    function residual(u)
        # Residual of the ODE at the nodes.
        du_dx = Dₓ * u                   # discrete u'
        d2u_dx2 = Dₓₓ * u                # discrete u''
        f = d2u_dx2 - ϕ.(x, u, du_dx)

        # Replace first and last values by boundary conditions.
        f[1] = g₁(u[1], du_dx[1]) / h
        f[n+1] = g₂(u[n+1], du_dx[n+1]) / h
        return f
    end

    u = levenberg(residual, init)
    return x, u[end]
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
  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.

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

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

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 ~/.julia/packages/FNCFunctions/Hk7fC/src/chapter04.jl:178

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¶

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

Preface

Collocation for linear problems

Preface

The Galerkin method