Ein Löser für 2-2-Bandmatrizen

The topic

Calculating spline curves requires the solution of linear equation systems, where the number of equations is close to the number of grid points. For cubic splines the matrices of the corresponding equation systems are occupied only along the diagonal and the subdiagonal on both sides. For such systems, simple solutions are common.

For higher order splinefunktionen or other border condition problems it may prove useful to be able to solve systems with two subdiagonals on each side (2-2-band matrices). Such a solution is derived and demonstrated below.

LU-Decomposition

The solution of the equation system is based upon the LU decomposition and the fact - the proof being omitted in this context - that the decomposition is again a 2-2 band matrix.

Let's consider an initial2-2 band matrix A and take a look at the first two steps of the decomposition:

Step 1

The first step transforms matrix A = A⁽¹⁾ into matrix A⁽²⁾

Which elements are modified and how?

• Row 1 remains unaltered
• Column 1 from row 2 results by dividing the elements a₂₁...a_n1 by a₁₁, and given the band structure, only a₂₁ and a₃₁ are considered, the others are zero
• The submatrix a⁽²⁾₂₂ bis a⁽²⁾₃₃ results by subtracting the dyadic vector product of vektors [a⁽²⁾₂₁, a⁽²⁾₃₁]^T mit [a⁽¹⁾₁₂, a⁽¹⁾₁₃] from [[a⁽¹⁾₂₂, a⁽¹⁾₂₃], [a⁽¹⁾₃₂, a⁽¹⁾₃₃]]
  a⁽²⁾₂₂ = a⁽¹⁾₂₂ – (a⁽¹⁾₁₂a⁽¹⁾₂₁)/a⁽¹⁾₁₁ usw.
• Column 4, coloured in orange in the image, also results from column 4 by subtraction of the product of the column 1 vector with element a14. The latter contains zero, leaving column 4 unaltered, and the same holds for the columns up to n.
• Analog for row 4, coloured in yellow, und the other rows down to n

In summary, only the elemente of sub matrix A⁽¹⁾_2..3,1..3 are modified.

Step 2

Step 2 transforms matrix A⁽²⁾ into A⁽³⁾ , and again only the elemetns in a 2x3 block are altered.

Which elements are modified and how?

• Row 2 remains unaltered
• Column 2, elements a⁽³⁾₃₂ and a⁽³⁾₄₂ result from a⁽²⁾₃₂ und a⁽²⁾₄₂ by dividing by a⁽²⁾₂₂ from step 1. Both manipulations need to be contained in the formula:
  a⁽³⁾₃₂ = (a⁽¹⁾₃₂ – (a⁽²⁾₁₂a⁽²⁾₃₁)/a⁽²⁾₁₁)/a⁽²⁾₂₂) usw.
• Sub matrix a⁽³⁾₃₃ to a⁽³⁾₄₄ results by subtracting the dyadic vektor produkts of vektors [a⁽³⁾₃₂, a⁽³⁾₄₂]^T and [a⁽²⁾₂₃, a⁽²⁾₂₄] from [[a⁽²⁾₃₃, a⁽²⁾₃₄], [a⁽²⁾₄₃, a⁽²⁾₄₄]]
  Element a⁽³⁾₃₃ needs special mention, obtained from a⁽²⁾₃₃. Element a⁽²⁾₃₃ cannot be stored intermediately, the calculation formula thus includes Step 1 of a⁽¹⁾₃₃ into a⁽²⁾₃₃. For the elements in row and column 4 the input values can be taken from the original matrix.

For the remaining steps k the input values are available either from A(k-1), that is in the calculation range, or from the originalmatrix A.

Example

Equation system Ax - b = 0

Matrix A  and condition vector b

The calculation schema implemented holds up to 20 equations, also less. No pivot elements are searched. In the input range of A columns a₂ and a₁ are the lower secondary diagonals, columns c₁ and c₂ the upper ones.

LU decomposition and vektor y

The LU decomposition of A in a format similar to the memory saving convolution of the full matrix form: The L-Matrix, designated M bezeichnet, has only ones in the diagonal which are not stored but contained in the calculation formulae.

Solution vector x

The solution of the equation system.