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, but 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 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 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}, \qquad \mathbf{Q} = \begin{bmatrix} q(x_0) & & \\ & \ddots & \\ & & q(x_{n}) \end{bmatrix},\\ \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

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

(10.4.9)

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.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
function [x, u] = bvplin(p, q, r, a, b, ua, ub,n)
% BVPLIN   Solve a linear boundary-value problem.
% Input:
%   p, q, r  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;

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.1 (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.10)

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

Julia

MATLAB

Python

Example 10.4.1

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

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

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

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 Demo 10.2.3, which exposed instability in the shooting method, in order to verify second-order convergence.

Example 10.4.2 (Convergence for a linear BVP)

The BVP is

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

(10.4.11)

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

Julia

MATLAB

Python

Example 10.4.2

λ = 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(size(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])
pretty_table(data, header = ["n", "inf-norm error"])

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

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
disp(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.2

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)

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}(h), \end{gather*}

(10.4.12)

where $\boldsymbol{\tau}$ is the truncation error of the finite differences (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 $h\to 0$ . 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¶

✍ For each boundary-value problem, verify that the given solution is correct. Then write out by hand for $n=3$ the matrices $\mathbf{D}_{xx}$ , $\mathbf{D}_x$ , $\mathbf{P}$ , and $\mathbf{Q}$ , and the vector $\mathbf{r}$ .
(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}$
⌨ For each of the cases in the previous exercise, use Function 10.4.1 to solve the problem with $n=60$ and make a plot of its error as a function of $x$ . Then, for each $n=10,20,40,\ldots,640$ , find the infinity norm of the error. Make a log-log plot of error versus $n$ and include a graphical comparison to second-order convergence.
⌨ Modify Function 10.4.1 to use spectral differentiation rather than second-order finite differences. For each of the cases in Exercise 1, solve the problem with $n=5,10,15,\ldots,40$ , finding the infinity norm of the error in each case. Make a log-linear plot of error versus $n$ .
⌨ Use Function 10.4.1 to solve Bessel’s equation,
$x^2 u'' + x u' + x^2 y = 0, \quad u(0.5)=1,\; u(8) = 0.$
(10.4.13)
Plot the solution for $n=100$ .
⌨ 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$ .
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.14)
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