
The first time I read the expression "orthogonal collocation" was in the article by Lee & Froment (2008) about the styrene reactor modeling. Turns out second-order differential equation systems appear in many Engineering topics, mostly related to transport phenomena. Thus, methods to solve them can be quite useful in practice. Once a process is mathematically described, it can be the object of numerical optimization problems, which, in this context, are usually related to process design and control.
If this content sparks your interest, take a look at the code repository of the Python framework used, also available on ResearchGate.
Orthogonal collocation was proposed by Villadsen and Stewart (1967) and can be very useful to solve systems of equations that look like this:

The math behind solutions can be quite complex, but, if you can formulate your system of equations like this, with appropriate boundary conditions, I believe OrthogonalCollocation (the Python framework) can do some magic.
Example 1
Let us get started with a simple second-order problem with only one dependent variable y.

Let us start by importing the main class of the framework, OrthogonalCollocation.
from collocation import OrthogonalCollocation
We must then create one function that returns zeros in internal points and another that returns zeros in the surface boundary (as the origin in x is assumed to be in symmetry conditions). Both can be either linear or nonlinear functions of x, y, y’, y”.
def fun_1(x, y, dy, d2y, k):
return d2y[0] + k * y[0] + 1
def bc_1(x, y, dy, d2y, k):
return dy[0] - 1
After defining the functions, the process is almost done, and we will solve the problem it in just three more lines of code (could have been just two). I’m using k=1 here, but we could have tried different numbers. Feel free to download the example notebook in my code repository and try different functions. Notice that, when performing collocations, scipy.optimize.root is used to solve systems of nonlinear equations implicit in the problems. All additional keywords in collocate are passed to root.
n_points = 6
# Create problem
problem_1 = OrthogonalCollocation(fun_1, bc_1, n_points, 1, x0=0.0, x1=1.0)
# Initial estimation
y01 = np.zeros([1, n_points + 1])
# Collocation using scipy.optimize.root in backend
problem_1.collocate(y01, args=k, method="hybr", tol=1e-6)
It is ready now! The collocation points are stored in a property y of our instance _problem1. Its shape is (m, n), so each line corresponds to a different variable and columns correspond to different collocation points. We can also calculate the interpolated polynomial values of y for other coordinates rather than just the collocation points. See for that I’ll use the interpolate method from our problem instance.
Now let us plot the results and compare them to the results obtained by the peer solver _scipy.integrate.solvebvp to assure the framework is working well.
fig, ax = plt.subplots(figsize=[6, 4], dpi=200)
x = np.linspace(0, 1, 50)
ax.plot(x, problem_1.interpolate(x)[0], color="indigo", alpha=0.5, linestyle=":")
ax.plot(problem_1.x, problem_1.y[0], color="indigo", marker=".", label="OC", linestyle=" ")
ax.plot(res_scipy1.x, res_scipy1.y[0], color="green", alpha=0.5, linestyle=":", label="scipy")
ax.set_ylabel("y")
ax.set_xlabel("x")
ax.legend()
fig.tight_layout()
plt.show()
And it looks like this:

Example 2
An example with a single dependent variable worked very well. In this section, let us try a system of equations with a nonlinear term ln.

Once again we formulate the equations:
def fun_2(x, y, dy, d2y, k1, k2):
return np.array([d2y[0] + k1 * y[1] + 1, d2y[1] + k2 * np.log(1 + y[0])])
def bc_2(x, y, dy, d2y, k1, k2):
return np.array([y[0], y[1] - 1])
I’m using k1=1 and k2=-1 here to make it more challenging. And the same second step…
problem_2 = OrthogonalCollocation(fun_2, bc_2, n_points, 1, x0=0.0, x1=1.0)
y02 = np.zeros([2, n_points + 1])
problem_2.collocate(y02, args=(k1, k2), method="hybr", tol=1e-6)
And now it looks like this:

And once again the solver returned great results!
In the following example, let us try a complex real-world problem.
The Styrene Reactor Problem
Remember I mentioned the Styrene reactor? So I have been working on this problem for some time and I have been able to develop an insightful study in which orthogonal collocation was used to model intraparticle diffusion. I believe this can give a good example of how a real-world application would look like.
Once fundamental equations were able to describe the process, the model could predict how much of the process inputs (reactants) were converted into desired outputs (styrene) for different conditions (such as temperature and pressure). Then it was the object of an optimization task "What are the most adequate operational conditions to convert inputs into desired outputs?".
If you are interested in details, I invite you to take a look at the article Simulation and optimization of axial-flow and radial-flow reactors for dehydrogenation of ethylbenzene into styrene based on a heterogeneous kinetic model.
In simple terms (as much as possible), there is a diffusional problem, and the concentrations of reactants inside the catalyst pellets are different from those on the surface, and so are effective reaction rates, thus we must calculate the effectiveness factors of reactions.
In my example notebook, I’ve tested different numbers of collocation points, described equations in detail, and compared results to scipy. The complete formulation is quite long for a short article, so I will present here what the main results look like for different numbers of collocation points.

Notice three collocation points produce results with already impressive precision, and results with six points are almost equal to nine. All these results took less than a quarter of computational time compared to the benchmark _scipy.integrate.solvebvp…
Conclusions
A simple yet robust framework for solving symmetric boundary value problems using orthogonal collocation was developed in Python. It is easy to use and was validated against peer solvers. A complex real-world problem was implemented, in which with very few collocation points results were remarkably accurate. The framework is available on this GIT repository.
References
Lee, W. J. & Froment, G. F., 2008. Ethylbenzene Dehydrogenation into Styrene: Kinetic Modeling and Reactor Simulation. Ind. Eng. Chem. Res., 47(23), pp. 9183–9194.
Leite, B., Costa, A. O. S. & Costa Junior, E. F., 2021. Simulation and optimization of axial-flow and radial-flow reactors for dehydrogenation of ethylbenzene into styrene based on a heterogeneous kinetic model. Chem. Eng. Sci., Volume 244, p. 116805.
Villadsen, J. & Stewart, W. E., 1967. Solution of boundary-value problems by orthogonal collocation. Chem. Eng. Sci., 22(11), pp. 1483–1501.