As I understand it, collocation method for partial differential equations is something akin to interpolation. First we characterize the solution space as a linear combination of some set of linearly independent functions $\phi_i(x)$. The appropriate choice of $\phi_i(x)$ depends on the problem. But assuming that the anticipated solution is sufficiently smooth, the simplest and easiest case to consider is the space of polynomials on $x$. Thus, the solution is of the form $u(x) = \sum \alpha_i \phi_i(x)$ where each $\phi_i(x)$ is a different (linearly independent) polynomial function.
The goal is to determine which coefficients $\alpha_i$ make the best approximation. We do so by substituting the linear combination directly into the partial differential equation. Once you do this, you obtain expressions in terms of the functions $\phi_i(x)$ (lets call this expression EQ1). Now, you need to force this expression to be exact at a finite number of points x. So, we pick a set of points $\{x_i\}$ and evaluate EQ1 at these points. This results in a linear system of equations whose unknowns are the $\alpha_i$'s.
To put the process into context, let's consider the 1 dimensional boundary value problem:
$U_{xx} = f$ on $\Omega = (0,1)$
$u=0$ on $\partial \Omega$
Suppose I expect the solution to be smooth enough that it can be characterized by polynomials. I can choose $\phi_i=x^i$. Then my linear combination would be $u(x)=\sum \alpha_i x^i$ for $i=0,...,n$. I can substitute this linear combination directly into the pde and obtain:
$$U_{xx} = f$$
$$\frac{\partial^2}{\partial x^2}[\sum_{i=1}^n \alpha_i \phi_i(x)]=f(x)$$
$$\sum_{i=0}^n \alpha_i \phi_i''(x)=f(x)$$ (lets call this last equation EQ1).
Now, we need to choose points in our domain $\Omega = (0,1)$ to force EQ1 to be true. Since there are n+1 coefficients ($\alpha_i$), we need to choose n+1 points. While theoretically we can choose any points in the domain, certain choices of points produce better accuracy than others. In my experience on smooth problems, the points are chosen as the n+1 gaussian quadrature points on (0,1). Keep in mind that this is probably not the optimal choice for all problems and that each specific problem may have a better choice of points. Once we know the points $x_k$, we evaluate EQ1 at each point $x_k$ and obtain:
$$\sum_{i=1}^n \alpha_i \phi_i''(x_k)=f(x_k)$$ for $k=0,...,n$
Thus, we have a system of equations for the unknowns $\alpha_i$, which we can solve by whatever method we see fit. Once we know the exact values of the $\alpha_i$'s, the final solution, again, is just their linear combination:
$u(x) = \sum_{i=0}^n\alpha_i \phi_i(x)$
Unfortunately I can't access the journal where your article is, so I can't say with certainty that this is how it works for an integro-differential equation. But this is the general framework of the collocation method.
Update
After considering the comments, I realized there are 2 other important conditions on the choice of functions $\phi_i(x)$:
- They must all satisfy the boundary conditions.
- Ideally, they should have small support.
My initial choice of basis functions $\phi_i(x)=x^i$ do not all satisfy the boundary condition at $u(1)=0$. Furthermore, polynomial functions are non-zero over a large region, leading to a dense matrix system which can be costly to solve. Often, cubic hermite functions are used in place of the usual polynomial functions because:
- They can naturally be adjusted to satisfy the boundary conditions.
- They lead to sparse linear systems of equations