Collocation Methods from Interpolation

Table of Contents

Intro #︎

Collocation methods are a way of approximately solving differential equations (both partial and ordinary). They are very simple to code and implement.

This article will cover collocation methods as an extension of simple interpolation. If you’re not familiar with interpolation theory, see Interpolation Theory 101. Even if you are, I’ll give a quick rundown of the parts of interpolation theory that you’ll need for this article:

Interpolation Fundamentals #︎

1. We want to find a function $g(x)$ such that it meets a set of $m+1$ constraints $g(x_i) = y_i$
2. We assume $g(x)$ to be of the form $g(x) = \sum_{i=0}^N c_i \phi_i(x)$ ¹
3. To ensure existence and uniqueness of the solution, we prescribe that the number of degrees-of-freedom (DOFs) $N+1$ equals the number of constraints $m+1$
4. Expanding $g(x)$ for each constraint results in a system of equations with $N+1$ DOFs and $m+1$ equations.

Normal interpolation sets constraints on the value of the function $g$ itself. I also introduced slope constraints, which sets a constraint on the slope of the polynomial rather than the polynomial itself.

Collocation is a further continuation of that; we set the constraint(s) that the function must satisfy a differential equation at select locations.

Collocation Overview #︎

In collocation, we’re doing the exact same procedure as in interpolation, but with different constraints. These constraints come in two forms:

1. Boundary conditions
2. Satisfying a differential equation

Boundary condition constraints look identical to normal interpolation constraints; we prescribe the value of our function (or of the function derivative) at some specific location.² The satisfaction of the differential equation is where things are a bit more interesting.

Collocation for a Generic Differential Equation #︎

For a given differential equation (DE), assume we can write it in the form:

$$\mathcal{L}(u) = f$$

for $\mathcal{L}$ the DE operator, $u$ the solution function, and $f$ some source function. For example, the Poisson PDE is given by $u’’(x) = f$, therefore $\mathcal{L}(u) = u’’$. Another, more complicated example, would be the Blasius Boundary layer equation, which is given as $2u’’’ + u’’ u = 0$. In this case, $\mathcal{L}(u) = 2u’’’ + u’’ u$.

The DE is subject to boundary conditions (and/or initial conditions). We will assume that the domain of the problem is $x = [0,1]$, and thus the boundaries are at $x=0$ and $x=1$. We assume two boundary conditions, though generally the number of boundary conditions for a general DE is determined by

$$ \begin{align*} u'(1) &= a \\ u''(1) &= b \end{align*} $$

Defining the Solution Function #︎

Following the steps from interpolation, we need to assume a form for our function. We’ll denote this function as $u^h$, which will be an approximate solution to the PDE. We assume the form:

$$u^h(x) = \sum_{n=0}^N c_n \phi_n(x)$$

where $c_n$ are our DOFs/coefficients and $\phi_n(x)$ are our basis functions. The choice of $\phi_n(x)$ is left arbitrary for this article. A few options are Legendre, Laguerre, and Chebyshev polynomials (which all have their own special properties, not discussed here), or non-polynomial basis functions, such as the Fourier series.

By defining $u=u^h$, we approximate the continuous PDE given above as: $$ \begin{align*} \mathcal{L}(u^h) &= f \\ u'^{\,h}(1) &= a \\ u''^{\,h}(1) &= b \end{align*} $$

Defining the Constraints #︎

Before we determine the number of DOFs ($N+1$), we need to know the number of constraints. I’ve already mentioned that the boundary conditions represent constraints, so each boundary condition represents a constraint. We’ll denote the number of boundary constraints as $m_b$.

But for “Satisfying the differential equation”, the PDE is defined for the entire domain, not discrete locations. This is equivalent to defining an infinite number of constraints, one for for every possible location within the domain. This is obviously not feasible.

Therefore, we choose a select number of locations in the domain to enforce the PDE. The locations that we use are called collocation points. We’ll denote the number of collocation points as $m_c$.

So the total number of constraints are the number of boundary conditions plus the number of collocation points. Since we need the number of DoFs to match the number of constraints (see Interpolation Theory 101), we get:

$$N + 1 = m_b + m_c $$

Determining the Collocation Points #︎

We still need to determine how many collocation points and where on the domain they should be. This is a more nuanced discussion, so I’ve moved a more detailed discussion of it to Appendix A.

For right now, we’ll assume the number of collocation points is arbitrary and that the points are evenly distributed across the domain. We denote the set of collocation point locations as:

$$\left\{\xi_j \right\}_{j=1}^{m_ c}$$

Creating the System of Equations #︎

To create the system of equations, just like with interpolation, we simply express each constraint in terms of $u^h$ and expand the definition of it. So for our boundary conditions, we end up with:

$$ \begin{align*} u'^{\,h}(1) = a \quad &\Rightarrow \quad c_0 \phi_0'(1) + c_1 \phi_1'(1) + \dots + c_N \phi_N'(1) = a \\ u''^{\,h}(1) = b \quad &\Rightarrow \quad c_0 \phi_0''(1) + c_1 \phi_1''(1) + \dots + c_N \phi_N''(1) = b \end{align*} $$

The slope constraints section of my interpolation theory post has a derivation for $u'^{\, h}(x) = \sum_{n=0}^N c_n \phi_n'(x)$.

Next, we do the same for the collocation points. This forms a set of $m_c$ equations:

$$ \begin{align*} \mathcal{L}(u^h(\xi_1)) &= f(\xi_1) \\ \mathcal{L}(u^h(\xi_2)) &= f(\xi_2) \\ \vdots \\ \mathcal{L}(u^h(\xi_{m_c})) &= f(\xi_{m_c}) \\ \end{align*} $$

Matrix Form #︎

If $\mathcal{L}$ is linear, we can express the PDE as:

$$\mathcal{L}(\sum_{n=0}^N c_n \phi_n(x)) = \sum_{n=0}^N c_n \mathcal{L}(\phi_n(x))$$

An example of a linear operator would be the Poisson equation $\mathcal{L}(u) = u’’$, since the derivative operator is a linear operator itself.

We can take the coefficients $c_n$ and move them into a separate vector $c$ and can form a $(N+1) \times (N+1)$:

$$ \begin{bmatrix} \phi_0(1) & \phi_1(1) & \cdots & \phi_N(1) \\[8pt] \phi'_0(1) & \phi'_1(1) & \cdots & \phi'_N(1) \\[8pt] \mathcal{L}(\phi_0(\xi_1)) & \mathcal{L}(\phi_1(\xi_1)) & \cdots & \mathcal{L}(\phi_N(\xi_1)) \\[8pt] \vdots & \vdots & \ddots & \vdots \\[8pt] \mathcal{L}(\phi_0(\xi_{m_c})) & \mathcal{L}(\phi_1(\xi_{m_c})) & \cdots & \mathcal{L}(\phi_N(\xi_{m_c})) \end{bmatrix} \begin{bmatrix} c_0 \\[8pt] c_1 \\[8pt] c_2 \\[8pt] \vdots \\[8pt] c_N \\[8pt] \end{bmatrix} = \begin{bmatrix} a \\[8pt] b \\[8pt] f(\xi_1) \\[8pt] \vdots \\[8pt] f(\xi_{m_c}) \\[8pt] \end{bmatrix} $$

Residual form #︎

If $\mathcal{L}$ is a nonlinear differential operator, we cannot bring the sum over coefficients out of the operator argument. An example of a nonlinear differential operator would be the Blasius equation $\mathcal{L}(u) = 2u’’’ + u’’ u$. The $u’’ u$ term is specifically what makes the differential operator non-linear.

Since we cannot separate the coefficients from a nonlinear operator argument, we cannot form a matrix system of equations. Instead, we need to create a residual function, which can then be used by a nonlinear equation solver to find our solution.

To obtain the residual function, we simply take the system of equations defined above and shift all the terms to the left-hand side such that they all equal zero. Since the coefficients $c_n$ completely define $u^h$, we often say that the residual is a function of the coefficients. Doing this, we obtain:

$$ R(u^h) = R \left( \begin{bmatrix} c_0 \\[8pt] c_1 \\[8pt] c_2 \\[8pt] \vdots \\[8pt] c_N \\[8pt] \end{bmatrix} \right) = \begin{bmatrix} \sum_{i=0}^N c_i \phi_i(1) - a \\[8pt] \sum_{i=0}^N c_i \phi'_i(1) -b \\[8pt] \mathcal{L}\left(u^h(\xi_1)\right) - f(\xi_1) \\[8pt] \vdots \\[8pt] \mathcal{L}\left(u^h(\xi_{m_c})\right) - f(\xi_{m_c}) \\[8pt] \end{bmatrix} $$

Example PDE: Poisson #︎

Let’s apply this to the 1-D Poisson equation for a domain of $x \in [0,1]$:

$$ \begin{align*} u''(x) &= f(x) \quad x \in [0,1] \\ u(0) &= 1 \\ u'(1) &= 2 \end{align*} $$

where the last two equations are our boundary conditions.

We take the same general form of $u^h$ as in the general case:

$$u^h(x) = \sum_{n=0}^N c_n \phi_n(x)$$

By defining $u=u^h$, we approximate the continuous PDE given above as:

$$ \begin{align*} u''^{\,h}(x) &= f(x) \quad x \in [0,1] \\ u^h(0) &= 1 \\ u'^{\,h}(1) &= 2 \end{align*} $$

Creating the System of Equations #︎

For our boundary conditions, we end up with:

$$ \begin{align*} u^h(0) = 1 \quad &\Rightarrow \quad c_0 \phi_0(0) + c_1 \phi_1(0) + \dots + c_N \phi_N(0) = 1 \\ u'^{\,h}(1) = 2 \quad &\Rightarrow \quad c_0 \phi_0'(1) + c_1 \phi_1'(1) + \dots + c_N \phi_N'(1) = 2 \end{align*} $$

For the PDE itself, we just evaluate the PDE approximated with $u^h$ at the collocation points:

$$ \begin{align*} u''^{\,h}(\xi_1) = \quad &\Rightarrow \quad c_0 \phi_0''(\xi_1) + c_1 \phi_1''(\xi_1) + \dots + c_N \phi_N''(\xi_1) = f(\xi_1) \\ u''^{\,h}(\xi_2) = \quad &\Rightarrow \quad c_0 \phi_0''(\xi_2) + c_1 \phi_1''(\xi_2) + \dots + c_N \phi_N''(\xi_2) = f(\xi_2) \\ \vdots \\ u''^{\,h}(\xi_{m_c}) = \quad &\Rightarrow \quad c_0 \phi_0''(\xi_{m_c}) + c_1 \phi_1''(\xi_{m_c}) + \dots + c_N \phi_N''(\xi_{m_c}) = f(\xi_{m_c}) \\ \end{align*} $$

Since the Poisson operator is linear, we can form the following matrix problem:

$$ \begin{bmatrix} \phi_0(0) & \phi_1(0) & \dots & \phi_N(0) \\ \phi_0'(1) & \phi_1'(1) & \dots & \phi_N'(1) \\ \phi_0''(\xi_1) & \phi_1''(\xi_1) & \dots & \phi_N''(\xi_1) \\ \phi_0''(\xi_2) & \phi_1''(\xi_2) & \dots & \phi_N''(\xi_2) \\ \vdots & \vdots & \ddots & \vdots \\[8pt] \phi_0''(\xi_{m_c}) & \phi_1''(\xi_{m_c}) & \dots & \phi_N''(\xi_{m_c}) \\ \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_N \\ \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \\ f(\xi_1) \\ f(\xi_2) \\ \vdots \\ f(\xi_{m_c}) \\ \end{bmatrix} $$

And tada! You have the collocation method for the Poisson PDE. All you need to do is solve this matrix equation for the DOFs $c_n$ and then reconstruct $u^h$ to have your solution!

Appendix A: Determining Collocation Points Continued #︎

Note we have yet to determine the number of collocation points or the location of the collocation points.

How many? #︎

Choosing the number of collocation points $m_c$ is completely arbitrary. Even if $u^h$ exactly satisfies the PDE at the collocation points, there are no guarantees how well (or poorly) $u^h$ satisfies the PDE for the rest of the domain. What we do know is that the more points you have, the more accurate your solution becomes (generally). How many points is “good enough” is highly problem dependent.

Beyond accuracy though, increasing the number of points also increases the size of your system of equations. Depending on your choice of collocation locations and basis functions, the system of equations you end up with may be infeasible to solve.

Where should the collocation points be located? #︎

This is also an arbitrary choice; they can be wherever you want them to be. But this choice can have significant implications on the stiffness of your resulting system of equations.

The general advice is to choose locations that are the roots of some $m_c -1$ order orthogonal polynomial. The reason for this choice is beyond the scope of the article, but it is intimately related to Gauss quadrature (which also uses the roots of orthogonal polynomials).

A good general choice to place collocation points at Chebyshev nodes, which are the roots of Chebyshev polynomials.