Collocation for linear problems - Fundamentals of Numerical Computation

Let us now devise a numerical method based on finite differences for the linear TPBVP

u'' + p(x)u' + q(x)u = r(x), \quad u(a)=\alpha,\; u(b)=\beta.

(10.4.1)

The first step is to select nodes $x_0=a < x_1< \cdots < x_n=b$ . For finite differences these will most likely be equally spaced, while for spectral differentiation, they will be Chebyshev points.

Rather than solving for a function, we will solve for a vector of its approximate values at the nodes:

\mathbf{u} = \begin{bmatrix} u_0 \\ u_1 \\ \vdots \\ u_{n-1} \\ u_n \end{bmatrix} \approx \begin{bmatrix} \hat{u}(x_0) \\ \hat{u}(x_1) \\ \vdots \\ \hat{u}(x_{n-1}) \\ \hat{u}(x_{n}) \end{bmatrix},

(10.4.2)

where $\hat{u}$ is the exact solution of (10.4.1). If we so desire, we can use interpolation to convert the values $(x_i,u_i)$ into a function after the discrete solution is found.

10.4.1Collocation¶

Having defined values at the nodes as our unknowns, we impose approximations to the ODE at the same nodes. This approach is known as collocation. Derivatives of the solution are found (approximately) using differentiation matrices. For example,

\begin{bmatrix} \hat{u}'(x_0) \\[1mm] \hat{u}'(x_1) \\ \vdots \\ \hat{u}'(x_n) \end{bmatrix} \approx \mathbf{u}' = \mathbf{D}_x \mathbf{u},

(10.4.3)

with an appropriately chosen differentiation matrix $\mathbf{D}_x$ . Similarly, we define

\begin{bmatrix} \hat{u}''(x_0) \\[1mm] \hat{u}''(x_1) \\ \vdots \\ \hat{u}''(x_n) \end{bmatrix} \approx \mathbf{u}'' = \mathbf{D}_{xx} \mathbf{u},

(10.4.4)

with $\mathbf{D}_{xx}$ chosen in accordance with the node set.

The discrete form of (10.4.1) at the $n+1$ chosen nodes is

\mathbf{u}'' + \mathbf{P}\mathbf{u}' + \mathbf{Q}\mathbf{u} = \mathbf{r},

(10.4.5)

where

\begin{split} \mathbf{P} = \begin{bmatrix} p(x_0) & & \\ & \ddots & \\ & & p(x_{n}) \end{bmatrix}, \quad \mathbf{Q} = \begin{bmatrix} q(x_0) & & \\ & \ddots & \\ & & q(x_{n}) \end{bmatrix}, \quad \mathbf{r} = \begin{bmatrix} r(x_0) \\ \vdots \\ r(x_n) \end{bmatrix}. \end{split}

(10.4.6)

If we apply the definitions of $\mathbf{u}'$ and $\mathbf{u}''$ and rearrange, we obtain

\mathbf{L} \mathbf{u} = \mathbf{r}, \qquad \mathbf{L} = \mathbf{D}_{xx} + \mathbf{P}\mathbf{D}_x + \mathbf{Q},

(10.4.7)

which is a linear system of $n+1$ equations in $n+1$ unknowns.

We have not yet incorporated the boundary conditions. Those take the form of the additional linear conditions $u_0=\alpha$ and $u_n=\beta$ . We might regard this situation as an overdetermined system, suitable for linear least-squares. However, it’s usually preferred to impose the boundary conditions and collocation conditions exactly, so we need to discard two of the collocation equations to keep the system square. The obvious candidates for deletion are the collocation conditions at the two endpoints. We may express these deletions by means of a matrix that is an $(n+1)\times(n+1)$ identity with the first and last rows deleted:

\mathbf{E} = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ \vdots & & \ddots & & \vdots \\ 0 & \cdots & 0 & 1 & 0 \end{bmatrix} = \begin{bmatrix} \mathbf{e}_1^T \\ \vdots \\ \mathbf{e}_{n-1}^T \end{bmatrix},

(10.4.8)

where, as always, $\mathbf{e}_k$ is the $k$ th column (here starting from $k=0$ ) of an identity matrix. The product $\mathbf{E} \mathbf{A}$ deletes the first and last rows of $\mathbf{A}$ , leaving a matrix that is $(n-1)\times(n+1)$ . Similarly, $\mathbf{E}\mathbf{r}$ deletes the first and last rows of $\mathbf{r}$ .

Finally, we note that $\hat{u}(a)= \mathbf{e}_0^T\mathbf{u}$ and $\hat{u}(b)= \mathbf{e}_n^T\mathbf{u}$ , so the linear system including both the ODE and the boundary condition collocations is

\underbrace{ \begin{bmatrix} \mathbf{e}_0^T \\[1mm] \mathbf{E}\mathbf{L} \\[1mm] \mathbf{e}_n^T \end{bmatrix} }_{\mathbf{A}} \mathbf{u} = \underbrace{ \begin{bmatrix} \alpha \\[1mm] \mathbf{E}\mathbf{r} \\[1mm] \beta \end{bmatrix} }_{\mathbf{b}}.

(10.4.9)

We have arrived at an $(n+1)\times (n+1)$ linear system $\mathbf{A}\mathbf{u}=\mathbf{b}$ , which we can solve for $\mathbf{u}$ using any of the methods we have investigated.

Example 10.4.1 (Collocation for a linear BVP)

Let’s look at the discretization of the linear BVP

u'' = u - x u', \quad u(0)=-2, \; u(2)=3,

(10.4.10)

on an equispaced grid with $n=4$ points. This gives $h=(2-0)/4 = 1/2$ .

Writing the ODE as $u'' + x u' - u = 0$ , we identify $p(x)=x$ , $q(x)=-1$ , and $r(x)=0$ . Hence,

\mathbf{P} = \begin{bmatrix} 0 & & & & \\ & 1/2 & & & \\ & & 1 & & \\ & & & 3/2 & \\ & & & & 2 \end{bmatrix}, \quad \mathbf{Q} = -\mathbf{I}_5, \quad \mathbf{r} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix}.

(10.4.11)

Using (10.3.6) and (10.3.10) as the differentiation matrices, we find

\begin{align*} \mathbf{L}& = \frac{1}{1/4} \begin{bmatrix} 2 & -5 & 4 & -1 & \\ 1 & -2 & 1 & & \\ & 1 & -2 & 1 & \\ & & 1 & -2 & 1 \\ & -1 & 4 & -5 & 2 \end{bmatrix} \\ & + \mathbf{P} \cdot \frac{1}{1/2} \begin{bmatrix} -3/2 & 2 & -1/2 & & \\ -1/2 & 0 & 1/2 & & \\ & -1/2 & 0 & 1/2 & \\ & & -1/2 & 0 & 1/2 \\ & & 1/2 & -2 & 3/2 \end{bmatrix} \\ & + \mathbf{Q} \\ & = \begin{bmatrix} 7 & -20 & 16 & -4 & 0 \\ 7/2 & -9 & 9/2 & 0 & 0 \\ 0 & 3 & -9 & 5 & 0 \\ 0 & 0 & 5/2 & -9 & 11/2 \\ 0 & -4 & 18 & -28 & 13 \end{bmatrix}. \end{align*}

(10.4.12)

The system $\mathbf{L} \mathbf{u} = \mathbf{r}$ represents the discretized ODE. When we modify the system to include the boundary conditions, we arrive at

\mathbf{A} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 7/2 & -9 & 9/2 & 0 & 0 \\ 0 & 3 & -9 & 5 & 0 \\ 0 & 0 & 5/2 & -9 & 11/2 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} -2 \\ 0 \\ 0 \\ 0 \\ 3 \end{bmatrix}.

(10.4.13)

The system $\mathbf{A} \mathbf{u} = \mathbf{b}$ is the collocated BVP.

10.4.2Implementation¶

Our implementation of linear collocation is Function 10.4.1. It uses second-order finite differences but makes no attempt to exploit the sparsity of the matrices. It would be trivial to change the function to use spectral differentiation.

Algorithm 10.4.1 (bvplin)

Julia

MATLAB

Python

Solution of a linear boundary-value problem

bvplin.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 linear boundary-value problem

bvplin.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
function [x, u] = bvplin(p, q, r, a, b, ua, ub,n)
% BVPLIN   Solve a linear boundary-value problem.
% Input:
%   p, q, r  define u'' + pu' + qu = r (functions)
%   a, b     endpoints of the domain (scalars)
%   ua       value of u(a)
%   ub       value of u(b)
%   n        number of subintervals (integer)
% Output:
%   x       collocation nodes (vector, length n+1)
%   u       solution at nodes (vector, length n+1)

    [x, Dx, Dxx] = diffmat2(n, [a, b]);

    P = diag(p(x));
    Q = diag(q(x));
    L = Dxx + P * Dx + Q;    % ODE expressed at the nodes
    r = r(x);    

    % Replace first and last rows using boundary conditions. 
    I = speye(n+1); 
    A = [ I(:, 1)'; L(2:n, :); I(:, n+1)' ];

    f = [ ua; r(2:n); ub ];

    % Solve the system.
    u = A \ f;
end

Solution of a linear boundary-value problem

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

    Return vectors of the nodes and the solution values.
    """
    x, Dx, Dxx = diffmat2(n, xspan)

    P = np.diag(p(x))
    Q = np.diag(q(x))
    L = Dxx + P @ Dx + Q  # ODE expressed at the nodes

    # Replace first and last rows using boundary conditions.
    I = np.eye(n + 1)
    A = np.vstack([I[0], L[1:-1], I[-1]])
    b = np.hstack([lval, r(x[1:-1]), rval])

    # Solve the system.
    u = np.linalg.solve(A, b)

    return x, u

Example 10.4.2 (Solving a linear BVP)

We solve the linear BVP

u'' - (\cos x) u' + (\sin x) u = 0, \quad u(0)=1, \; u\left(\frac{3\pi}{2}\right)=\frac{1}{e}.

(10.4.14)

Its exact solution is $u(x) = e^{\sin x}$ .

Julia

MATLAB

Python

Example 10.4.2

exact = x -> exp(sin(x));

The problem is presented above in our standard form, so we can identify the coefficient functions in the ODE. Each should be coded as a function.

p = x -> -cos(x);
q = sin;
r = x -> 0;      # function, not value

We solve the BVP and compare the result to the exact solution.

x, u = FNC.bvplin(p, q, r, [0, 3π / 2], 1, exp(-1), 30);

using Plots
plot(exact, 0, 3π / 2, layout = (2, 1), label = "exact")
scatter!(x, u, m = :o,
    subplot=1,  label="numerical",
    yaxis=("solution"),
    title="Solution of a linear BVP")
plot!(x, exact.(x) - u, subplot = 2, xaxis = L"x", yaxis = ("error"))

Example 10.4.2

exact = @(x) exp(sin(x));

The problem is presented above in our standard form, so we can identify the coefficient functions in the ODE. Each should be coded as a function.

p = @(x) -cos(x);
q = @(x) sin(x);
r = @(x) 0*x;      % not a scalar

We solve the BVP and compare the result to the exact solution.

[x, u] = bvplin(p, q, r, 0, 3*pi/2, 1, exp(-1), 30);

clf,  subplot(2, 1, 1)
plot(x, u)
ylabel('solution')
title('Solution of a linear BVP')
subplot(2, 1, 2)
plot(x, exact(x) - u, 'o-')
ylabel('error')

Example 10.4.2

exact = lambda x: exp( sin(x) )

The problem is presented above in our standard form, so we can identify the coefficient functions in the ODE. Each should be coded as a function.

p = lambda x: -cos(x)
q = sin
r = lambda x: 0 * x    # must be a function

We solve the BVP and compare the result to the exact solution.

x, u = FNC.bvplin(p, q, r, [0, pi/2], 1, exp(1), 25)

subplot(2, 1, 1)
plot(x, u)
ylabel("solution"),  title("Solution of the BVP")

subplot(2, 1, 2)
plot(x, exact(x) - u, "-o")
ylabel("error");

10.4.3Accuracy and stability¶

We revisit Example 10.2.3, which exposed instability in the shooting method, in order to verify second-order convergence.

Example 10.4.3 (Convergence for a linear BVP)

The BVP is

u'' - \lambda^2 u = \lambda^2, \quad u(0)=-1, \; u(1)=0,

(10.4.15)

with exact solution $\sinh(\lambda x)/\sinh(\lambda) - 1$ .

Julia

MATLAB

Python

Example 10.4.3

λ = 10
exact = x -> sinh(λ * x) / sinh(λ) - 1;

The following functions define the ODE.

p = x -> 0
q = x -> -λ^2
r = x -> λ^2;

We compare the computed solution to the exact one for increasing $n$ .

n = 5 * [round(Int, 10^d) for d in 0:0.25:3]
err = zeros(length(n))
for (k, n) in enumerate(n)
    x, u = FNC.bvplin(p, q, r, [0, 1], -1, 0, n)
    err[k] = norm(exact.(x) - u, Inf)
end
data = (n = n[1:4:end], err = err[1:4:end])
@pt :header = ["n", "inf-norm error"] data

Each factor of 10 in $n$ reduces error by a factor of 100, which is indicative of second-order convergence.

plot(n, err, m = :o,
    label = "observed",
    xaxis = (:log10, L"n"),
    yaxis = (:log10, "inf-norm error"),
    title = "Convergence for a linear BVP")
plot!(n, 0.25 * n .^ (-2), l = (:dash, :gray), label = "2nd order")

Example 10.4.3

lambda = 10;
exact = @(x) sinh(lambda * x) / sinh(lambda) - 1;

The following functions define the ODE.

p = @(x) zeros(size(x));            
q = @(x) -lambda^2 * ones(size(x));
r = @(x) lambda^2 * ones(size(x));

We compare the computed solution to the exact one for increasing $n$ .

p = @(x) zeros(size(x));            
q = @(x) -lambda^2 * ones(size(x));
r = @(x) lambda^2 * ones(size(x));
n = 2 * round(10.^(1:0.25:3)');
err = zeros(size(n));
for k = 1:length(n)
    [x, u] = bvplin(p, q, r, 0, 1, -1, 0, n(k));
    err(k) = norm(exact(x) - u, Inf);
end
table(n, err, variableNames = ["n", "inf-norm error"])

Each factor of 10 in $n$ reduces error by a factor of 100, which is indicative of second-order convergence.

clf,  loglog(n, err, 'o-')
hold on, loglog(n, n.^(-2), 'k--')
xlabel('n'),  ylabel('max error')
title('Convergence for a linear BVP') 
legend('obs. error', '2nd order')

Example 10.4.3

lamb = 10
exact = lambda x: sinh(lamb * x) / sinh(lamb) - 1

The following functions define the ODE.

p = lambda x: zeros(size(x))
q = lambda x: -(lamb**2) * ones(len(x))
r = lambda x: lamb**2 * ones(len(x))

We compare the computed solution to the exact one for increasing $n$ .

N = array([int(2 * 10**d) for d in arange(1, 3.1, 0.25)])
err = zeros(len(N))
results = PrettyTable(["n", "error"])
for k, n in enumerate(N):
    x, u = FNC.bvplin(p, q, r, [0, 1], -1, 0, n)
    err[k] = norm(exact(x) - u, inf)
    results.add_row([n, err[k]])
print(results)

+------+------------------------+
|  n   |         error          |
+------+------------------------+
|  20  | 0.0037471006137431706  |
|  35  | 0.0012307958808570607  |
|  63  | 0.0003848512001197735  |
| 112  | 0.00012208705546823406 |
| 200  | 3.8311986798711395e-05 |
| 355  | 1.216083231025955e-05  |
| 632  | 3.8374952795461326e-06 |
| 1124 | 1.2132645536011566e-06 |
| 2000 | 3.8320567119498605e-07 |
+------+------------------------+

Each factor of 10 in $n$ reduces error by a factor of 100, which is indicative of second-order convergence.

loglog(N, err, "-o", label="observed")
loglog(N, 1 / N**2, "--", label="2nd order")
xlabel("$n$"),  ylabel("max error")
legend(),  title("Convergence of finite differences");

If we write the solution $\mathbf{u}$ of Equation (10.4.9) as the exact solution minus an error vector $\mathbf{e}$ , i.e., $\mathbf{u} = \hat{\mathbf{u}} - \mathbf{e}$ , we obtain

\begin{gather*} \mathbf{A} \hat{\mathbf{u}} - \mathbf{A} \mathbf{e} = \mathbf{b}, \\ \mathbf{e} = \mathbf{A}^{-1} \left[ \mathbf{A} \hat{\mathbf{u}} - \mathbf{b} \right] = \mathbf{A}^{-1} \boldsymbol{\tau}(n), \end{gather*}

(10.4.16)

where $\boldsymbol{\tau}$ is the truncation error of the derivative discretizations (except at the boundary rows, where it is zero). It follows that $\|\mathbf{e}\|$ vanishes at the same rate as the truncation error if $\| \mathbf{A}^{-1}\|$ is bounded above as $n\to \infty$ . In the present context, this property is known as stability. Proving stability is too technical to walk through here, but stability is guaranteed under some reasonable conditions on the BVP.

10.4.4Exercises¶

Exercise 10.4.1

✍ For each boundary-value problem, verify that the given solution is correct. Using $n=3$ and the differentiation matrices in (10.3.6) and (10.3.10), write out $\mathbf{L}$ and $\mathbf{r}$ from (10.4.7) and $\mathbf{A}$ and $\mathbf{b}$ from (10.4.9).

(a) $u'' + u = 0, \quad u(0) =0, \; u(3) = \sin 3$

Solution: $u(x) = \sin x$

(b) $u'' - \frac{3}{x} u' + \frac{4}{x^2} u = 0, \quad u(1) =0,\; u(4) = 32 \log 2$

Solution: $u(x) = x^2 \log x$

(c) $u'' - \left(x+\frac{1}{2}\right)^{-1}\, u' + 2\left(x+\frac{1}{2}\right)^{-2}\, u = 10\left(x+\frac{1}{2}\right)^{-4}, \quad u\left(x+\frac{1}{2}\right)=1,\; u\left(x+\frac{5}{2}\right) = \frac{1}{9}$

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

Exercise 10.4.5

⌨ The Airy equation is $u''=x u$ . Its solution is exponential for $x>0$ and oscillatory for $x<0$ . The exact solution is given by $u=c_1 \operatorname{Ai}(x) + c_2 \operatorname{Bi}(x)$ , where Ai and Bi are Airy functions. In Julia they are computed by airyai and airybi, respectively.

(a) Suppose that $u(-10) =-1$ , $u(2) =1$ . By setting up and solving a $2\times 2$ linear system, find numerical values for $c_1$ and $c_2$ . Plot the resulting exact solution.

(b) Use Function 10.4.1 with $n=120$ to find the solution with the boundary conditions in part (a). In a 2-by-1 subplot array, plot the finite-difference solution and its error. (The solution is not very accurate.)

(c) Repeat part (b) with $n=800$ .

Exercise 10.4.6

Consider the boundary-value problem $\epsilon u''+(1+\epsilon)u'+u =0$ over $x\in (0,1)$ , with $u(0)=0$ , $u(1)=1$ . As the parameter ε is decreased, the solution gets a thin region of high activity near $x=0$ called a boundary layer.

(a) ✍ Verify that the exact solution to the problem is

u(x) = \frac{e^{-x}-e^{-x/\epsilon}}{e^{\,-1}-e^{\,-1/\epsilon}}.

(10.4.18)

On one graph, plot $u(x)$ for $\epsilon=\frac{1}{4},\frac{1}{16},\frac{1}{64}$ .

(b) ⌨ Define $N(\epsilon)$ as the smallest integer value of $n$ needed to make the max-norm error of the result of Function 10.4.1 less than 10^-4. For each of the values $\epsilon = \frac{1}{2},\frac{1}{4},\frac{1}{8}\ldots,\frac{1}{64}$ , estimate $N(\epsilon)$ by starting with $n=50$ and incrementing by 25 until the measured error is sufficiently small.

(c) ⌨ Plot the error as a function of $x$ for $\epsilon=\frac{1}{64}$ and $n=N(\epsilon)$ . Compare the peak of the error to the graph from part (a).

(d) ⌨ Develop a hypothesis for the leading-order behavior of $N(\epsilon)$ . Plot the observed $N(\epsilon)$ and your hypothesis together on a log-log plot.

(e) ✍ Finite-difference errors depend on the solution as well as on $n$ . Given that this error decreases as $O(n^{-2})$ , what does your hypothesis for $N(\epsilon)$ suggest about the behavior of the error for fixed $n$ as $\epsilon\to 0$ ?

Preface

Differentiation matrices

Preface

Nonlinearity and boundary conditions