matrix equations for a spline

matrix equation for an open curve

The matrix equation to be solved is $$ \begin{pmatrix} b_{0} & c_0 & & & 0 \\ a_{1} & b_{1} & c_1 & & \\ & a_{2} & b_{2} & \ddots & \\ & & \ddots & \ddots & c_{n-1} \\ 0 & & & a_n & b_n \end{pmatrix} \begin{pmatrix} P_{1,0} \\ P_{1,1} \\ P_{1,2} \\ \vdots \\ P_{1,n} \end{pmatrix} = \begin{pmatrix} r_{0} \\ r_{1} \\ r_{2} \\ \vdots \\ r_{n} \end{pmatrix} \quad\quad\quad (12) $$ in which $$ \begin{align} b_{0} &= 2 &\\ c_{0} &= w_0/w_1 &\\ a_{i} &= w_i^2 & i \in \{1,\dots, n-1\} \\ b_{i} &= 2w_{i-1}(w_{i-1}+w_i) & i \in \{1,\dots, n-1\} \\ c_{i} &= (w_{i-1})^2w_i/w_{i+1} & i \in \{1,\dots, n-1\} \\ a_{n} &= 1 &\\ b_{n} &= 2 w_{n-1}/w_n& \end{align} $$ and $$ \begin{align} r_0 &= K_{0} + K_1 \left( 1+ {w_{0} \over w_{1}} \right) \\ r_i &= K_{i}\left( w_{i-1} + w_{i}\right)^2 + K_{i+1} w_{i-1}^2 \left( 1+ {w_{i} \over w_{i+1}} \right) & i \in \{1,\dots, n-1\} \\ r_n &= K_n \left( 1+2 {w_{n-1} \over w_{n}} \right) \end{align} $$ (corrected on 2018-08-16 after a comment by R.B.A. van Hoek)
The equation must be solved 2 times: 1 time for the \(x\) coordinates, and 1 time for the \(y\) coordinates of \(\mathbf P_{1,i}\).

matrix equation for a closed curve

In this case, \(n\) knots will be labeled \(0,1,\dots,n-1\). $$ \begin{pmatrix} b_{0} & c_0 & & & a_0 \\ a_{1} & b_{1} & c_1 & & \\ & a_{2} & b_{2} & \ddots & \\ & & \ddots & \ddots & c_{n-2} \\ c_{n-1} & & & a_{n-1} & b_{n-1} \end{pmatrix} \begin{pmatrix} P_{1,0} \\ P_{1,1} \\ P_{1,2} \\ \vdots \\ P_{1,n-1} \end{pmatrix} = \begin{pmatrix} r_{0} \\ r_{1} \\ r_{2} \\ \vdots \\ r_{n-1} \end{pmatrix} \quad\quad\quad (13)$$ Again we have $$ \begin{align} a_{i} &= w_i^2 & i \in \{0,\dots, n-1\} \\ b_{i} &= 2w_{i-1}(w_{i-1}+w_i) & i \in \{1,\dots, n-1\} \\ c_{i} &= (w_{i-1})^2w_i/w_{i+1} & i \in \{1,\dots, n-2\} \\ r_i &= K_{i}\left( w_{i-1} + w_{i}\right)^2 + K_{i+1} w_{i-1}^2 \left( 1+ {w_{i} \over w_{i+1}} \right) & i \in \{1,\dots, n-2\} \\ \end{align} $$ Of course, we use index \(0\) in stead of \(n\), so the formulas yield: $$ \begin{align} b_{0} &= 2w_{n-1}(w_{n-1}+w_0) \\ c_{0} &= (w_{n-1})^2w_0/w_{1} \\ c_{n-1} &= (w_{n-2})^2w_{n-1}/w_{0} \\ r_0 &= K_{0}\left( w_{n-1} + w_{0}\right)^2 + K_{1} w_{n-1}^2 \left( 1+ {w_{0} \over w_{1}} \right) \\ r_{n-1} &= K_{n-1}\left( w_{n-2} + w_{n-1}\right)^2 + K_{0} w_{n-2}^2 \left( 1+ {w_{n-1} \over w_{0}} \right) \end{align} $$

Gaussian elimination

The equation (12) can be solved using Gaussian elimination. Because the matrix for the open curve is tri-diagonal the Thomas algorithm can be used.
For a looped curve (equation (13)), I will describe the algorithm below.

The non-zero entries of the matrix will be stored in arrays a, b and c with indices 0 till n-1.
The right-hand-side is the array r with indices 0 till n-1.
During Gaussian elimination, the elements in the most right column of the matrix (under \(a_0\)) will change and I will store them in the array u with indices 0 till n-2.
The elements on the lowest row will also change, but become zero again. So, in stead of an array L with indices 0 till n-2, I use a single variable L.
In this way, the algorithm becomes:

u[0] = a[0]
L = c[n-1]
for (i = 0; i < n-3; i++)
{
	m1 = a[i+1]/b[i];
	b[i+1] -= m1 * c[i];
	// last column 
	u[i+1] = -m1 * u[i]

	// last row : 
	m2 = L/b[i]  
	b[n-1] -= m2 * u[i]
	L = - m2 * c[i] 

	// right-hand side
	r[i+1] -= m1 * r[i];
	r[n-1] -= m2 * r[i]
}
After this, the matrix equation has become $$ \begin{pmatrix} b_{0} & c_0 & & & & u_0 \\ 0 & b_{1} & c_1 & & & u_1\\ & \ddots & \ddots & \ddots & & \vdots\\ & & 0 & b_{n-3} & c_{n-3} & u_{n-3}\\ & & & a_{n-2} & b_{n-2} & c_{n-2} \\ 0 & & & L & a_{n-1} & b_{n-1} \end{pmatrix} \begin{pmatrix} P_{1,0} \\ P_{1,1} \\ P_{1,2} \\ \vdots \\ P_{1,n-1} \end{pmatrix} = \begin{pmatrix} r_{0} \\ r_{1} \\ r_{2} \\ \vdots \\ r_{n-1} \end{pmatrix} $$ So the next steps must be (note that i = n-3) 
m1 = a[i+1]/b[i];
b[i+1] -= m1 * c[i];
// last column 
c[i+1] -= m1 * u [i]

// last row 
m2 = L/b[i]  
b[n-1] -= m2 * u[i]
a[n-1] -= m2 * c[i]

// right-hand side
r[i+1] -= m1 * r[i];
r[n-1] -= m2 * r[i]
which yields $$ \begin{pmatrix} b_{0} & c_0 & & & & u_0 \\ 0 & b_{1} & c_1 & & & u_1\\ & \ddots & \ddots & \ddots & & \vdots\\ & & 0 & b_{n-3} & c_{n-3} & u_{n-3}\\ & & & 0 & b_{n-2} & c_{n-2} \\ 0 & & & 0 & a_{n-1} & b_{n-1} \end{pmatrix} \begin{pmatrix} P_{1,0} \\ P_{1,1} \\ P_{1,2} \\ \vdots \\ P_{1,n-1} \end{pmatrix} = \begin{pmatrix} r_{0} \\ r_{1} \\ r_{2} \\ \vdots \\ r_{n-1} \end{pmatrix} $$ After 
i = n-2
m = a[i+1]/b[i];
b[i+1] -= m * c[i];
r[i+1] -= m * r[i];
this becomes $$ \begin{pmatrix} b_{0} & c_0 & & & & u_0 \\ 0 & b_{1} & c_1 & & & u_1\\ & \ddots & \ddots & \ddots & & \vdots\\ & & 0 & b_{n-3} & c_{n-3} & u_{n-3}\\ & & & 0 & b_{n-2} & c_{n-2} \\ 0 & & & & 0 & b_{n-1} \end{pmatrix} \begin{pmatrix} P_{1,0} \\ P_{1,1} \\ P_{1,2} \\ \vdots \\ P_{1,n-1} \end{pmatrix} = \begin{pmatrix} r_{0} \\ r_{1} \\ r_{2} \\ \vdots \\ r_{n-1} \end{pmatrix} $$ Then the backward substitution becomes easy: 
x[n-1] = r[n-1]/b[n-1];
x[n-2] = (r[n-2] - c[n-2] * x[n-1]) / b[n-2];
for (i = n - 3; i >= 0; --i)
{	x[i] = (r[i] - c[i] * x[i+1]-u[i]*x[n-1]) / b[i];
}
Previous