From Function Spaces to Matrix Representation

Integration by parts applied to $\int e^{2x}(3x^2+x-1)\,dx$ is a familiar routine: differentiate the polynomial, integrate the exponential, repeat until the polynomial runs out. The routine always ends with an answer of the same shape, $e^{2x}$ times a quadratic. Chapter 5 of Rida Abu-Sokon's book turns that observation into linear algebra. If the functions involved form a finite-dimensional space that differentiation maps into itself, then differentiation is a matrix, and the integral can be found by solving one linear system. This explainer builds the idea from the vector-space vocabulary up.

Functions as vectors

A set of functions is a vector space when sums and scalar multiples of its members stay in the set. Polynomials of degree at most two form one: adding $1+x$ and $x^2$ gives another polynomial of that kind. The functions $a\sin\omega x+b\cos\omega x$ form another. A linear combination of functions $\phi_1,\dots,\phi_m$ is any expression $c_1\phi_1+\dots+c_m\phi_m$; their span is the set of all such combinations; and they form a basis of the space they span when none of them is a combination of the others. The number of basis functions is the dimension. The space $\operatorname{span}\{1,x,x^2\}$ has dimension three even though it contains infinitely many functions.

Once a basis is fixed, every function in the space is recorded by its list of coefficients. Collect the basis into a column $\Phi(x)=(\phi_1(x),\dots,\phi_m(x))^{T}$; then $f(x)=u^{T}\Phi(x)$, where $u\in\mathbb{R}^m$ is the coordinate vector of $f$. For example, $3x^2+x-1$ has coordinates $u=(-1,1,3)^{T}$ in the basis $(1,x,x^2)$.

Closure under differentiation

A space $V$ is closed under differentiation if $f\in V$ implies $f'\in V$. The space $\operatorname{span}\{1,x,x^2\}$ is closed, and so is $\operatorname{span}\{\sin\omega x,\cos\omega x\}$. The space $\operatorname{span}\{x\}$ is not, because $(x)'=1$ is not a multiple of $x$; neither is $\operatorname{span}\{1/x\}$. Closure is the property that makes everything below possible, because it lets differentiation act inside a fixed, finite list of coordinates.

The differentiation matrix

If $V$ is closed, the derivative of each basis function is a combination of basis functions, $\phi_i'=\sum_j\Omega_{ij}\phi_j$. In vector form,

So on coordinates, differentiation is the matrix $\Omega^{T}$. For the two bases above (the book lists these and others as "ready matrices"):

Check the first one on $3x^2+x-1$: $\Omega^{T}(-1,1,3)^{T}=(1,6,0)^{T}$, which is the coordinate vector of $6x+1$. The trigonometric matrix is $\omega$ times a quarter-turn rotation: differentiating $\cos$ and $\sin$ rotates the coordinate vector by $90^\circ$ and scales it by $\omega$.

Why polynomials give nilpotent matrices

Differentiation lowers the degree of a polynomial by one. On polynomials of degree at most $m-1$, applying it $m$ times gives zero for every input, so the matrix satisfies $(\Omega^{T})^m=0$. A matrix with a vanishing power is called nilpotent; its only eigenvalue is $0$, and it is strictly triangular in the monomial basis. This is not a quirk of the basis: any finite-dimensional space of polynomials closed under differentiation forces nilpotency, because a nonzero eigenvalue $\lambda$ would need a polynomial with $p'=\lambda p$, and only $e^{\lambda x}$ satisfies that.

The spectrum of $\Omega$ therefore records what kind of functions live in the space. Polynomials give eigenvalue $0$; $\cos\omega x,\sin\omega x$ give $\pm i\omega$; $e^{\alpha x}\cos\omega x,e^{\alpha x}\sin\omega x$ give $\alpha\pm i\omega$; $e^{\alpha x}\cosh\omega x,e^{\alpha x}\sinh\omega x$ give $\alpha\pm\omega$.

From an integral to a linear solve

Now take $f=u^{T}\Phi\in V$ and look for an antiderivative of $e^{bx}f(x)$ of the same shape, $F(x)=e^{bx}\Phi(x)^{T}z$ with an unknown constant vector $z$. The product rule and $\Phi'=\Omega\Phi$ give $F'(x)=e^{bx}\Phi(x)^{T}(bI+\Omega^{T})z$. Matching this with $e^{bx}\Phi(x)^{T}u$ and using linear independence of the basis gives a single linear system:

This boundary formula is what the book calls the Matrix Boundary Method (MBM). Its recipe has five steps: choose a basis and write $f=u^{T}\Phi$; build $\Omega$; invert $bI+\Omega^{T}$ once; form $F(x)=e^{bx}\Phi^{T}M(b)u$; evaluate $F(a)-F(0)$. Because $M(b)$ depends only on the basis and on $b$, one inversion serves every $f$ in the space.

For polynomials, nilpotency also explains the shape of $M(b)$. Since $(\Omega^{T}/b)^m=0$, the geometric series for the inverse terminates:

The powers $1/b,1/b^2,1/b^3$ are the successive steps of integration by parts, now produced by three terms of a finite series.

Because $M(2)$ does not depend on $f$, any other quadratic costs only a matrix–vector product. For $f=x^2$, $u=(0,0,1)^{T}$ gives $z=(\tfrac14,-\tfrac12,\tfrac12)^{T}$, so $\int_0^1e^{2x}x^2\,dx=\big[e^{2x}(\tfrac14-\tfrac12x+\tfrac12x^2)\big]_0^1=(e^2-1)/4\approx1.59726$. By hand, each new polynomial would mean a fresh round of integration by parts.

The trigonometric basis works the same way. With $b=1$ and $\omega=2$, $M=\frac15\begin{pmatrix}1&-2\\2&1\end{pmatrix}$; for $f=\sin 2x$, $u=(0,1)^{T}$ and $z=(-\tfrac25,\tfrac15)^{T}$, so $\int_0^\pi e^{x}\sin 2x\,dx=\big[e^{x}(-\tfrac25\cos 2x+\tfrac15\sin 2x)\big]_0^\pi=\tfrac25(1-e^{\pi})\approx-8.85628$.

The link to the Laplace transform

The truncated integral is called a truncated Laplace integral for a reason. Put $b=-s$ and let $a\to\infty$. If $s$ exceeds the real parts of all eigenvalues of $\Omega$, the boundary term at $a$ vanishes and only the value at $0$ survives:

On a finite-dimensional space closed under differentiation, the Laplace transform is therefore the resolvent $(sI-\Omega^{T})^{-1}$ of the differentiation matrix, read off at $x=0$. For polynomials the first row of $(sI-\Omega^{T})^{-1}$ is $(1/s,\,1/s^2,\,2/s^3)$, reproducing $\mathcal{L}\{1\}$, $\mathcal{L}\{x\}$ and $\mathcal{L}\{x^2\}$; for $(\cos\omega x,\sin\omega x)$ it gives $s/(s^2+\omega^2)$ and $\omega/(s^2+\omega^2)$. This is the same structure that appears in linear systems theory, where $(sI-A)^{-1}$ is the Laplace transform of the matrix exponential $e^{xA}$ that solves $y'=Ay$.

What the method is, and what it is not

The underlying facts are standard: coordinate representations of linear maps, the matrix of differentiation, and the method of undetermined coefficients for $y'+by=f$. What the book adds is a systematic packaging for truncated integrals $\int_0^a e^{bx}f(x)\,dx$, which are finite-interval versions of the Laplace integral, together with a catalogue of ready matrices and extensions to products of bases and to linear ordinary differential equations in later sections of the chapter. The method is exact, but it applies only to inputs that lie in a known finite-dimensional space closed under differentiation; it does not evaluate $\int_0^a e^{bx}\ln x\,dx$, for example.