Sign in Join
Wiki copytemporary

NumLib Documentation

Introduction

NumLib is a collection of routines for numerical mathematics. It was written by C.J.J.M. van Ginneken, W.J.P.M. Kortsmit, and L. Reij at the Technical University of Eindhoven (Netherlands) and donated to the Free Pascal project. Its beginnings date back to 1965. Originally written in Algol 60, ported to poorly formatted Pascal with strange procedure names, it is hard stuff, in particular because an official documentation is not readily available (well... -- Dutch documentation files in TeX can be found at http://www.stack.nl/~marcov/numlib/).

But after all, this library contains a series of very useful routines grouped into several units:

Therefore, it is the intention of this wiki article to create some useful documentation. Each routine described here was tested and its description will be accompanied by a working code example.

Note: This page heavily makes use of mathematical formulas written in LaTeX. In order to see these formulas in the browser JavaScript must be enabled.

Naming convention of routines

Would you imagine that a procedure named spears calculates the arcsin? The original NumLib routines limit the procedure/function names to 6 characters (maybe a limitation of Algol?). The first three characters always are identical with the name of the unit in which they are implemented. This leaves only three characters to differentiate between the routines...

Routines later added to the original library are given more conventional names, such as normaldist for calculation of the normal distribution.

Data types and declarations (unit typ)

Basic floating point type, ArbFloat

In NumLib, a new type ArbFloat is defined which is used as floating point type throughout the entire library. This allows to switch between IEEE Double (64bit) and Extended(80bit) (though big endian state is unknown at present).

At the time of this writing the $DEFINE ArbExtended is enabled, and therefore ArbFloat is identical to the standard type extended. In order to change the floating point type define or undefine ArbExtended and adapt the machineconstants to the type selected.

Basic integer type, ArbInt

Because in plain FPC mode the type Integer is 16-bits (for TP compatibility), and in Delphi 32-bits, all integers were changed to ArbInt. The basic idea is the same as with ArbFloat.

Vectors

NumLib often passes vectors to procedures as one ArbFloat plus some integer value. This ArbFloat value is the first value of an array containing the vector, the integer is the count of elements used from this array. Internally an array type is mapped over the first value in order to access the other array elements. This way both 0- and 1-based arrays as well as static and dynamic arrays can be passed to the functions and procedures of NumLib.

pascal
procedure DoSomething(var AVector: ArbFloat; n: Integer);
type
  TLocalVector = array[0..MaxElements] of ArbFloat;
var
  i: Integer;
  x: TLocalVector absolute AVector;
begin
  for i:=0 to n-1 do
    // do something with x[i]
end;

var
  x: array[1..100];   // Note: This array begins with index 1!
...
  DoSomething(x[1], 100);

Matrices

A matrix is a two-dimensional arrangement of numbers in m rows and n columns.

${A} = \left[ \begin{array}{cccc} a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\ a_{21} & a_{22} & a_{23} & \ldots & a_{2n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & a_{m3} & \ldots & a_{mn} \end{array} \right]$

The matrix elements a<sub>ij</sub> are addressed by row and column indexes, i and j, respectively. In mathematics, the index values usually start with 1 (i = 1..m, j = 1..n) while in computing the first index is preferred to be 0 (i = 0...m-1, j = 0...n-1).

Generic matrices

In this type of matrix, no assumption is made on symmetry and shape of the data arrangement.

A constant (real) matrix can be declared in Pascal as a two-dimensional array of floating point numbers: A: array[1..m, 1..n] of ArbFloat (or A: array[0..m-1, 0..n-1] of ArbFloat in 0-based index notation). Similarly, a constant complex matrix can be declared as an array[1..m, 1..n] of complex.

When a NumLib procedure needs a matrix the first (top left) matrix element must be passed as a var parameter, along with the number of rows and columns. It is even possible to allocate larger arrays as needed, and in this case the allocated row length (or number of columns) must be specified.

Alternatively a matrix also can be specified as a one-dimensional array. This works because all data are found within the same memory block. The array index, k must be converted to the 2D matrix indexes, i, j:

pascal
  k := i * n + j      <==> i = k div n,           j := k mod n          {for 0-based array and matrix indexes}
  k := (i-1) * n + j  <==> i := (k-1) div n + 1,  j := (k-1) mod n + 1  {for 1-based array and matrix indexes}

Using a one-dimensional array is also the correct way to dynamically allocate memory for a matrix. The first idea of using an array of array of ArbFloat will not work because this does not allocate a single memory block. Instead, individual memory blocks will be allocated for the inner array (rows), and these pointers will be stored in the outer array - this way of storing 2D data is not compatible with the requirements of NumLib.

Here is a simple example how to create, populate and read a matrix dynamically

math
A = \left[
  \begin{array}{ccc}
  1 & 2 & 3 \\
  4 & 5 & 6
  \end{array}
  \right]
pascal
var
  A: array of Arbfloat;  // Note: the index in a dynamic array is 0-based.
  ij: Integer;
  i, j: Integer;
  p: ^ArbFloat;
  m, n: Integer;
begin
  m := 2;
  n := 3;
  SetLength(A, m, n);
  A[0] := 1;  A[1] := 2;  A[2] := 3;  
  A[3] := 4;  A[4] := 5;  A[5] := 6;

  for i := 0 to m-1 do begin
    for j := 0 to n-1 do begin
      ij := i * n + j;
      Write(A[ij]);
    end;
    WriteLn;
  end;

  // or
  p := @A[0];   // pointer to first array element
  for i:=0 to m - 1 do begin
    for j := 0 to n - 1 do begin
      WriteLn(p^);
      inc(p);   // advance pointer to next array element
    end;
    WriteLn;
  end;

Symmetric matrices

The elements of a symmetric matrix are mirrored at the main diagonal, i.e. a<sub>ij</sub> = a<sub>ji</sub>. In other words: a symmetric matrix is equal to its transpose: A = A<sup>T</sup>.

NumLib contains some routine optimized for such matrices. They require the matrix being stored like a generic matrix in a rectangular m x n array or a linear m*n array. However, only the half of the matrix below the diagonal must be populated (but it can be, of course).

Symmetric positive definite matrices

A matrix A is positive definite if the product d<sup>T</sup> A d is a positive number for any non-zero vector d. A matrix is symmetric positive definite if, additionally, A<sup>T</sup> = A.

This type of matrices has some interesting properties (determinant is never 0, eigenvalues are real and positive). Therefore NumLib offers specialized routines for many matrix problems with symmetric positive matrices which are superior to those for general matrices.

The matrix is stored in the conventional way like a rectangular 2D array of m rows and n columns or a linear 1D array with m**n* elements.

Band matrices

A n x n band matrix has zero elements everywhere except for a band along the main diagonal characterized by a left and right bandwidth (l and r, respectively).

The NumLib procedures for band matrices assume that the matrix stores only the elements of the band. It is expected the the matrix is passed as a 1D array. Starting at the top/left corner point of the matrix (element A<sub>11</sub>) the diagonal elements are collected by running along the rows of the matrix, from left to right and from top to bottom. In total, there will be n (l+1+r) - (r(r+1) + l(l+1)) div 2 elements. The calling procedure must allocate a 1D array of at least this size for the diagonal elements. Note that a band matrix cannot be passed as a 2D array.

As an example, suppose this 7 x 7 matrix with left bandwidth 2 and right bandwidth 1:

$A = \begin{bmatrix} A_{11} & A_{12} & 0 & 0 & 0 & 0 & 0 \\ A_{21} & A_{22} & A_{23} & 0 & 0 & 0 & 0 \\ A_{31} & A_{32} & A_{33} & A_{34} & 0 & 0 & 0 \\ 0 & A_{42} & A_{43} & A_{44} & A_{45} & 0 & 0 \\ 0 & 0 & A_{53} & A_{54} & A_{55} & A_{56} & 0 \\ 0 & 0 & 0 & A_{64} & A_{65} & A_{66} & A_{67} \\ 0 & 0 & 0 & 0 & A_{75} & A_{76} & A_{77} \end{bmatrix}$

The 1D data array to be passed to the NumLib routines consists of these elements

pascal
const 
  arr: array[1..24] of ArbFloat = (    
    A11, A22,
    A21, A22, A23,
    A31, A32, A33, A34,
         A42, A43, A44, A45,
              A53, A54, A55, A56,
                   A64, A65, A66, A67,
                        A75, A76, A77
  );

Symmetric band matrices

In a symmetric matrix, the matrix elements can be mirrored at the main diagonal of the matrix, i.e. A<sub>ij</sub> = A<sub>ji</sub>. All elements out side a band of width w to the left and right of the main diagonal are zero.

Like with a generic band matrix, NumLib assumes that the only the elements of the diagonals are passed as a 1D array to the procedures dedicated to this matrix type. Since in a symmetric band matrix the left and right band widths are equal it is assumed that only the left band and the main diagonal are stored. Again the relevant diagonal elements are collected by running across the rows of the matrix from left to right and from top to bottom, stopping in each row at the main diagonal. In total, there will be n (w+1) - w(w+1) div 2 elements in the data array.

As an example consider this 7 x 7 symmetric band matrix with band width 2:

$A = \begin{bmatrix} A_{11} & A_{21} & A_{31} & 0 & 0 & 0 & 0 \\ A_{21} & A_{22} & A_{32} & A_{42} & 0 & 0 & 0 \\ A_{31} & A_{32} & A_{33} & A_{43} & A_{53} & 0 & 0 \\ 0 & A_{42} & A_{43} & A_{44} & A_{54} & A_{64} & 0 \\ 0 & 0 & A_{53} & A_{54} & A_{55} & A_{65} & A_{75} \\ 0 & 0 & 0 & A_{64} & A_{65} & A_{66} & A_{76} \\ 0 & 0 & 0 & 0 & A_{75} & A_{76} & A_{77} \end{bmatrix}$

The 1D array for the NumLib procedures is

pascal
const
  arr: array[1..18] of ArbFloat = (
    A11,
    A21, A22,
    A31, A32, A33,
         A42, A43, A44,
              A53, A54, A55,
                   A64, A65, A66,
                        A75, A76, A77
  );

Tridiagonal matrices

A tridiagonal matrix has zero elements everywhere except for the main diagonal and the diagonal immediately above and below it. This way the matrix can be stored more efficiently than a generic matrix because all the empty elements can be dropped. The elements of an n x n square tridiagonal matrix can be stored as three vectors

$A= \begin{bmatrix} d[0] & u[0] & & \dots & 0 \\ l[1] & d[1] & u[1] & & \\ & l[2] & d[2] & u[2] & \\ \vdots & & \ddots & \ddots & \ddots \\ & & l[n-1] & d[n-1] & u[n-1] \\ 0 & & & l[n] & d[n] \end{bmatrix}$

Symmetric tridiagonal matrices

In these tridiagonal matrices, the subdiagonal and superdiagonal elements are equal. Consequently, NumLib requests only two, instead of three, one-dimensional arrays:

$A= \begin{bmatrix} d[0] & c[1] & & \dots & 0 \\ c[1] & d[1] & c[2] & & \\ & c[2] & d[2] & c[3] & \\ \vdots & & \ddots & \ddots & \ddots \\ & & c[n-1] & d[n-1] & c[n] \\ 0 & & & c[n] & d[n] \end{bmatrix}$

Complex numbers

NumLib defines the type complex as a Turbo Pascal object. Real and imaginary parts are denoted as xreal and imag, respectively. Methods for basic operations are provided:

pascal
type
  complex  = object
    xreal, imag : ArbFloat;
    procedure Init(r, i: ArbFloat);    // Initializes the number with real and imaginary parts r and i, resp
    procedure Add(c: complex);         // Adds another complex number 
    procedure Sub(c: complex);         // Subtracts a complex number
    function  Inp(z:complex):ArbFloat; // Calculates the inner product (a.Re * b.Re + a.Im * b.Im)
    procedure Conjugate;               // Conjugates the value (i.e. negates the imaginary part)
    procedure Scale(s: ArbFloat);      // Multiplies real and imaginary part by the real value s
    function  Norm: ArbFloat;          // Calculates the "norm", i.e. sqr(Re) + sqr(Im) 
    function  Size: ArbFloat;          // Size = sqrt(Norm) 
    function  Re: ArbFloat;            // Returns the real part of the complex number
    procedure Unary;                   // Negates real and imaginary parts
    function  Im: ArbFloat;            // Returns the imaginary part of the complex number
    function  Arg: ArbFloat;           // Calculates the phase phi of the complex number z = x + y i = r exp(i phi)
    procedure MinC(c: complex);        // Converts itself to the minimum of real and imaginary parts of itself and another value, c
    procedure MaxC(c: complex);        // Converts itself to the maximum of real and imaginary ports of itself and another value, c
    procedure TransF(var t: complex);
  end;

The real and imaginary parts of a complex number can also be retrieved by these stand-alone functions:

pascal
  function Re(c: complex): ArbFloat;
  function Im(c: complex): ArbFloat;

Basic operations with matrices and vectors (unit omv)

Inner product of two vectors

The result of the inner product of two vectors a = [a<sub>1</sub>, a<sub>2</sub>, ... a<sub>n</sub>] and b = [b<sub>1</sub>, b<sub>2</sub>, ... b<sub>n</sub>] (also called "dot product" or "scalar product") is a scalar, i.e. a single number. It is calculated as the sum of the element-by-element products of the vectors. Both vectors must have the same lengths.

$\mathbf{a} \cdot \mathbf{b} = \sum_{i=1}^n a_i b_i = a_1 b_1 + a_2 b_2 + \dots + a_n b_n$

NumLib provides the function omvinp for calculation of the inner product:

pascal
function omvinp(var a, b: ArbFloat; n: ArbInt): ArbFloat;

Example

Calculate the dot product of the vectors a = [0, 1, 2, 2, -1] and b = [3, -1, -2, 2, -1].

pascal
program dot_product;

uses
  typ, omv;

const
  n = 5;
  a: array[0..n-1] of ArbFloat = (0, 1, 2, 2, -1);
  b: array[0..n-1] of Arbfloat = (3, -1, -2, 2, -1);

var
  ab: ArbFloat;
  i: Integer;

begin
  ab := omvinp(a[0], b[0], n);

  Write('Vector a = [');
  for i:= Low(a) to High(a) do
    Write(a[i]:5:0);
  WriteLn('  ]');

  Write('Vector b = [');
  for i:= Low(b) to High(b) do
    Write(b[i]:5:0);
  WriteLn('  ]');

  Write('Dot product a b = ');
  WriteLn(ab:0:0);
end.

Product of two matrices

The product of a m x n matrix A and a n x p matrix B is a m x p matrix C in which the elements are calculated as the inner product of the row vectors of A and the column vectors of B (see here for more details):

$C_{ij} = \sum_{k=0}^n A_{ik} B_{kj}$

The NumLib function to perform this multipication is omvmmm ("mmm" = multiply matrix by matrix):

pascal
procedure omvmmm(var a: ArbFloat; m, n, rwa: ArbInt;
                 var b: ArbFloat; p, rwb: ArbInt;
                 var c: ArbFloat; rwc: ArbInt);

Example

Calculate the product of the matrices

$A = \begin{bmatrix} 3 & \ 5 & 4 & 1 & -4 \\ -2 & \ 3 & 4 & 1 & 0 \\ 0 & \ 1 & -1 & -2 & 5 \end{bmatrix} \ \ \ \ \ B = \begin{bmatrix} 3 & 4 & -2 & 0 \\ 0 & -2 & 2 & 1 \\ -1 & 0 & 3 & 4 \\ -2 & -3 & 0 & 1 \\ 1 & 3 & 0 & -1 \end{bmatrix}$

pascal
program matrix_multiplication;

uses
  typ, omv;

const
  m = 3;
  n = 5;
  p = 4;

  A: array[1..m, 1..n] of ArbFloat = (
    ( 3,  5,  4,  1, -4),
    (-2,  3,  4,  1,  0),
    ( 0,  1, -1, -2,  5)
  );
  B: array[1..n, 1..p] of ArbFloat = (
    ( 3,  4, -2,  0),
    ( 0, -2,  2,  1),
    (-1,  0,  3,  4),
    (-2, -3,  0,  1),
    ( 1,  3,  0, -1)
  );

var
  C: array[1..m, 1..p] of ArbFloat;
  i, j: Integer;

begin
  // Calculate product
  omvmmm(
    A[1,1], m, n, n,
    B[1,1], p, p,
    C[1,1], p
  );

  // Print A
  WriteLn('A = ');
  for i:= 1 to m do begin
    for j := 1 to n do
      Write(A[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;

  // Print B
  WriteLn('B = ');
  for i:= 1 to n do begin
    for j := 1 to p do
      Write(B[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;

  // Print C
  WriteLn('C = A B = ');
  for i:= 1 to m do begin
    for j := 1 to p do
      Write(C[i, j]:10:0);
    WriteLn;
  end;
end.

Product of a matrix with a vector

A matrix A of size m x n multiplied with a vector b of length n yields another vector of length m. The elements of the result vector c are the dot products of the matrix rows with b. Note that the length of vector b must be equal to the number of columns of matrix A.

$`\begin{align*} \mathbf{c} = A\ \mathbf{b}= \left[ \begin{array}{cccc} a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\ a_{21} & a_{22} & a_{23} & \ldots & a_{2n}\ \vdots & \vdots & \vdots & \ddots & \vdots\ a_{m1} & a_{m2} & a_{m3} & \ldots & a_{mn} \end{array} \right] \ \ \left[ \begin{array}{c} b_1\ b_2\ b_3\ \vdots\ b_n \end{array} \right]

\left[ \begin{array}{c} a_{11}b_1+a_{12}b_2 + a_{13}b_3 + \cdots + a_{1n} b_n\ a_{21}b_1+a_{22}b_2 + a_{23}b_3 + \cdots + a_{2n} b_n\ \vdots\ a_{m1}b_1+a_{m2}b_2 + a_{m3}b_3 + \cdots + a_{mn} b_n\ \end{array} \right]. \end{align*}`$

NumLib's procedure for this task is called omvmmv ("mmv" = multiply a matrix with a vector):

pascal
procedure omvmmv(var a: ArbFloat; m, n, rwidth: ArbInt; var b, c: ArbFloat);

Example

Calculate the product

$\begin{align*} \mathbf{c} = A\ \mathbf{b}= \left[ \begin{array}{cccc} 3 & 5 & 4 & 1 & -4\\ -2 & 3 & 4 & 1 & 0\\ 0 & 1 & -1 & -2 & 5 \end{array} \right] \ \ \left[ \begin{array}{c} 3\\ 0\\ -1\\ -2\\ 1 \end{array} \right] \end{align*}$

pascal
program matrix_vector_multiplication;

uses
  typ, omv;

const
  m = 3;
  n = 5;
  p = 4;

  A: array[1..m, 1..n] of ArbFloat = (
    ( 3,  5,  4,  1, -4),
    (-2,  3,  4,  1,  0),
    ( 0,  1, -1, -2,  5)
  );
  b: array[1..n] of ArbFloat = (
    3, 0, -1, -2, 1
  );

var
  c: array[1..m] of ArbFloat;
  i, j: Integer;

begin
  // Calculate product c = A b
  omvmmv(
    A[1,1], m, n, n,
    b[1],
    c[1]
  );

  // Print A
  WriteLn('A = ');
  for i:= 1 to m do begin
    for j := 1 to n do
      Write(A[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;

  // Print vector b
  WriteLn('b = ');
  for i:= 1 to n do
    Write(b[i]:10:0);
  WriteLn;
  WriteLn;

  // Print result vector c
  WriteLn('c = A b = ');
  for i:= 1 to m do
    Write(c[i]:10:0);
  WriteLn;
end.

Transpose matrix

The transpose matrix A<sup>T</sup> of matrix A is obtained by flipping rows and columns:

$\begin{align*} {A} = \left[ \begin{array}{cccc} a_{11} & a_{12} & a_{13} & \ldots & a_{1n}\\ a_{21} & a_{22} & a_{23} & \ldots & a_{2n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{m1} & a_{m2} & a_{m3} & \ldots & a_{mn} \end{array} \right] \ \ \ {A^T} = \left[ \begin{array}{cccc} a_{11} & a_{21} & \ldots & a_{m1}\\ a_{12} & a_{22} & \ldots & a_{m2}\\ a_{13} & a_{23} & \ldots & a_{m3} \\ \vdots & \vdots & \ddots & \vdots \\ a_{1n} & a_{2n} & \ldots & a{mn} \end{array} \right] \end{align*}$

Use the procedure omvtrm to perform this operation with NumLib:

pascal
procedure omvtrm(
  var a: ArbFloat; m, n, rwa: ArbInt; 
  var c: ArbFloat; rwc: ArbInt
);

Example

Obtain the transpose of the 2x4 matrix

${A} = \left[ \begin{array}{cccc} 1 & 2 & 3 & 4 \\ 11 & 22 & 33 & 44 \end{array} \right]$

pascal
program transpose;

uses
  typ, omv;

const
  m = 2;
  n = 4;

  A: array[1..m, 1..n] of ArbFloat = (
    ( 1,  2,  3,  4),
    (11, 22, 33, 44)
  );

var
  C: array[1..n, 1..m] of ArbFloat;
  i, j: Integer;

begin
  // Transpose A
  omvtrm(
    A[1,1], m, n, n,
    C[1,1], m
  );

  // Print A
  WriteLn('A = ');
  for i:= 1 to m do begin
    for j := 1 to n do
      Write(A[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;

  // Print C
  WriteLn('AT = ');
  for i:= 1 to n do begin
    for j := 1 to m do
      Write(C[i, j]:10:0);
    WriteLn;
  end;
end.

Norm of vectors and matrices

A norm is a function which assigns a positive "length" to a vector a or a matrix M.

In case of vectors, NumLib supports these norms:

$\|a\|_1 = \sum_{i=1}^n |{a_i}|$

$\|a\|_2 = \sqrt{\sum_{i=1}^n {a_i}^2}$

$\|a\|_\infty = \max({a_1}, {a_2}, ... {a_n})$

The norms of matrices, as implemented by NumLib, are usually defined by row or column vectors:

$\|M\|_1 = \max_{1 \le j \le {n}} \sum_{i=1}^m|M_{ij}|$

$\|M\|_\infty = \max_{1 \le i \le\ m} \sum_{j=1}^n |M_{ij}|$

$\|M\|_F = \sqrt{\sum_{i=1}^m \sum_{j=1}^n {M_{ij}}^2}$

These are the NumLib procedure for calculation of norms:

pascal
function omvn1v(var a: ArbFloat; n: ArbInt): ArbFloat;               // 1-norm of a vector
function omvn1m(var a: ArbFloat; m, n, rwidth: ArbInt): ArbFloat;    // 1-norm of a matrix
function omvn2v(var a: ArbFloat; n: ArbInt): ArbFloat;               // 2-norm of a vector
function omvnfm(Var a: ArbFloat; m, n, rwidth: ArbInt): ArbFloat;    // Frobenius norm of a matrix
function omvnmv(var a: ArbFloat; n: ArbInt): ArbFloat;               // Maximum infinite norm of a vector
function omvnmm(var a: ArbFloat; m, n, rwidth: ArbInt): ArbFloat;    // Maximum infinite norm of a matrix

Example

pascal
program norm;

uses
  typ, omv;

const
  n = 2;
  m = 4;
  a: array[0..m-1] of ArbFloat = (0, 1, 2, -3);
  b: array[0..m-1, 0..n-1] of Arbfloat = ((3, -1), (-2, 2), (0, -1), (2, 1));

var
  i, j: Integer;

begin
  Write('Vector a = [');
  for i:= Low(a) to High(a) do
    Write(a[i]:5:0);
  WriteLn('  ]');
  WriteLn('        1-norm of vector a: ', omvn1v(a[0], m):0:3);
  WriteLn('        2-norm of vector a: ', omvn2v(a[0], m):0:3);
  WriteLn('  maximum norm of vector a: ', omvnmv(a[0], m):0:3);
  WriteLn;

  WriteLn('Matrix b = ');
  for i:= 0 to m-1 do begin
    for j := 0 to n-1 do
      Write(b[i,j]:5:0);
    WriteLn;
  end;
  WriteLn('        1-norm of matrix b: ', omvn1m(b[0, 0], m, n, n):0:3);
  WriteLn('Forbenius norm of matrix b: ', omvnfm(b[0, 0], m, n, n):0:3);
  WriteLn('  maximum norm of matrix b: ', omvnmm(b[0, 0], m, n, n):0:3);

  ReadLn;
end.

Determinant of a matrix (unit det)

Unit det provides a series of routines for calculation of the determinant of a square matrix. Besides a procedure for general matrices there are also routines optimized for special matrix types (tridiagonal, symmetric).

In order to avoid possible overflows and underflows a factor may be split off during calculation, and the determinant is returned as a floating point number f times a multiplier 8<sup>k</sup>. In many cases, however, k = 0, and the determinant simply is equal to f.

Determinant of a standard matrix

The term "standard matrix" means here a matrix with is stored in the standard NumLib storage as described here, i.e. the matrix is stored as a 2D or 1D array of real numbers. The 2D array has at least n rows and n columns, the 1D array is dimensioned to contain at least n<sup>2</sup> elements.

Without making any assumptions on the properties of the matrix the determinant can be calculated by the procedure detgen. More efficient algorithms can be used, however, if the matrix has these particular properties:

pascal
procedure detgen(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt);    // --> generic matrix
procedure detgsy(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt);    // --> symmetric matrix
procedure detgpd(n, rwidth: ArbInt; var a, f: ArbFloat; var k, term: ArbInt);    // --> symmetric positive definite matrix

Example

Calculate the determinant of the 4x4 matrix

$A= \left[ \begin{array}{rrrr} 4 & 2 & 4 & 1 \\ 30 & 20 & 45 & 12 \\ 20 & 15 & 36 & 10 \\ 35 & 28 & 70 & 20 \end{array} \right]$

The matrix is positive definite because all elements are positive and the product b<sup>T</sup> A b cannot become zero or negative for any non-zero vector b. While the routine detgpd is the best for this matrix type the code below contains also commented-out calls to the other routines.

pascal
program determinant_spd_matrix;   // "spd" = symmetric positive definite

uses
  math, typ, det;

const
  n = 4;

const
  a: array[1..n, 1..n] of ArbFloat = (
    ( 4,  2,  4,  1)
    (30, 20, 45, 12)
    (20, 15, 36, 10)
    (35, 28, 70, 20)
  );

var
  f: ArbFloat;
  k: Integer;
  term: Integer = 0;
  i, j: Integer;

begin
  // write input matrix
  WriteLn('A = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(a[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;

  // Calculate determinant of symmetric positive definite matrix
  detgpd(n, n, a[1, 1], f, k, term);

  // or uncomment to use one of the other methods
  // detgen(n, n, a[1, 1], f, k, term);
  // detgsy(n, n, a[1, 1], f, k, term);

  if term = 1 then
    WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
  else
    WriteLn('Incorrect input');
end.

Determinant of a band matrix

If a matrix has the shape of a band matrix then its determinant is best calculated by using the procedure detgba

pascal
procedure detgba(n, l, r: ArbInt; var a, f: ArbFloat; var k, term:ArbInt);

Example

Calculate the determinant of the 5x5 band matrix with left bandwidth l = 2 and right bandwidth r = 1

math
A = \begin{bmatrix}
    1 & 2 & 0 & 0 & 0 \\
    2 & 1 & 1 & 0 & 0 \\
    1 & 3 & 0 & 1 & 0 \\
    0 & 1 & 1 & 1 & 1 \\
    0 & 0 & 1 & 1 & 1 
  \end{bmatrix}
pascal
program determinant_band_matrix;

uses
  math, typ, det;

const
  n = 5;
  L = 2;
  R = 1;
  nelements = n*(L+1+R) - (R*(R+1) + L*(L+1)) div 2;

const
  // Diagonal elements arranged in rows, elements outside of the band omitted!
  a: array[1..nelements] of Arbfloat = (
     1, 2,
     2, 1, 1,
     1, 3, 0, 1,
     1, 1, 1, 1,
     1, 1, 1
    );
var
  f: ArbFloat;
  term: Integer = 0;
  i, i1, i2, j, k: Integer;

begin
  // Write diagonal elements
  WriteLn('n = ', n);
  Writeln('Left band width L = ', L);
  WriteLn('Right band width R = ', R);
  Write('Diagonal elements of A = ', a[1]:0:0);
  for k := 2 to nelements do
    Write(a[k]:5:0);
  WriteLn;
  WriteLn;

  // write reconstructed band input matrix (not needed for calculation)
  WriteLn('Reconstructed A = ');
  k := 1;
  for j := 1 to n do begin
    i1 := j - L;
    i2 := i1 + L + R;
    for i:=1 to i1-1 do
      Write(0.0:5:0);
    for i := i1 to i2 do
      if (i >= 1) and (i <= n) then begin
        Write(a[k]:5:0);
        inc(k);
      end;
    for i:= i2 + 1 to n do
      if i <= n then
        Write(0.0:5:0);
    WriteLn;
  end;
  WriteLn;

  // Calculate determinant
  detgba(n, l, r, a[1], f, k, term);
  if term = 1 then
    WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
  else
    WriteLn('Incorrect input');
end.

Determinant of a symmetric positive definite band matrix

If a matrix has the shape of a band matrix and is symmetric positive definite then detgpb is the dedicated routine for calculation of its determinant.

pascal
procedure detgpb(n, w: ArbInt; var a, f: ArbFloat; var k, term:ArbInt);

Example

Calculate the determinant of the 7x7 symmetric band matrix with bandwidth w = 2 $A= \begin{bmatrix} 5 & -4 & 1 & & & & \\ -4 & 6 & -4 & 1 & & & \\ 1 & -4 & 6 & -4 & 1 & & \\ & 1 & -4 & 6 & -4 & 1 & \\ & & 1 & -4 & 6 & -4 & 1 \\ & & & 1 & -4 & 6 & -4 \\ & & & & 1 & -4 & 5 \end{bmatrix}$

pascal
program determinant_spd_band_matrix;     // "spd" = symmetric positive definite

uses
  math, typ, det;

const
  n = 7;
  w = 2;
  nelements = n * (w + 1) - (w * (w + 1)) div 2;

const
  // Band elements arragned in rows, elements outside of the band omitted!
  a: array[1..nelements] of Arbfloat = (
             5,
         -4, 6,
      1, -4, 6,
      1, -4, 6,
      1, -4, 6,
      1, -4, 6,
      1, -4, 5
    );

  function MatrixIndexToArrayIndex(i, j, n, w: Integer): Integer;

    function DiagElementIndex(i: Integer): Integer;
    var
      k: Integer;
    begin
      Result := 1;
      if i = 1 then
        exit;
      // Rows truncated at left
      for k := 2 to w do begin
        Result := Result + k;
        if k = i then
          exit;
      end;
      // full rows and those truncated at right
      for k := w+1 to n do begin
        Result := Result + w + 1;
        if k = i then
          exit;
      end;
      Result := n;
    end;

  var
    d: Integer;
  begin
    if j > i then begin
      Result := MatrixIndexToArrayIndex(j, i, n, w);
      exit;
    end;

    Result := DiagElementIndex(i);
    if (i = j) then
      exit;

    d := i - j;
    if d > w then
      Result := -1
    else begin
      dec(Result, d);
      if (Result < 1) then
        Result := -1
    end;
  end;

var
  i, j, k: Integer;
  f: ArbFloat;
  term: Integer = 0;
begin
  // Write diagonal elements
  WriteLn('n = ', n);
  Writeln('(One-sided) band width W = ', w);
  Write('Diagonal elements of A = ', a[1]:0:0);
  for k := 2 to nelements do
    Write(a[k]:3:0);
  WriteLn;
  WriteLn;

  // write reconstructed band input matrix (not needed for calculation)
  WriteLn('Reconstructed A = ');
  k := 1;
  for i := 1 to n do begin
    for j := 1 to n do begin
      k := MatrixIndexToArrayIndex(i, j, n, w);
      if k = -1 then
        Write(0.0:3:0)
      else
        Write(a[k]:3:0);
    end;
    WriteLn;
  end;
  WriteLn;

  // Calculate determinant
  detgpb(n, w, a[1], f, k, term);
  if term = 1 then
    WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
  else
    WriteLn('Incorrect input');
end.

Determinant of a tridiagonal matrix

Calculation of the determinant of a tridiagonal matrix is done by the procedure detgtr:

pascal
procedure detgtr(n: ArbInt; var l, d, u, f: ArbFloat; var k, term:ArbInt);

Example Calculate the determinant of the matrix

$A= \left( \begin{array}{rrrr} 1 & 1 & & \\ 2 & 1 & 1 & \\ & 1 & 1 & 0 \\ & & 1 & 2 \end{array} \right)$

pascal
program determinant_tridiagonal_matrix;

{$mode objfpc}{$H+}

uses
  math, typ, det;

const
  n = 4;

const
  u: array[0..n-2] of ArbFloat = (-1, -1, -1   );       // superdiagonal
  d: array[0..n-1] of ArbFloat = ( 2,  2,  2,  2);      // main diagonal
  l: array[1..n-1] of ArbFloat = (    -1, -1, -1);      // subdiagonal

var
  f: ArbFloat;
  k: Integer;
  term: Integer = 0;
  i, j: Integer;

begin
  // write input matrix
  WriteLn('A = ');
  for i := 0 to n-1 do begin
    for j := 0 to i-2 do
      Write(0.0:5:0);
    if i > 0 then
      Write(l[i]:5:0);
    Write(d[i]:5:0);
    if i < n-1 then
      Write(u[i]:5:0);
    for j:=i+2 to n-1 do
      Write(0.0:5:0);
    WriteLn;
  end;
  WriteLn;

  // Calculate determinant
  detgtr(n, l[1], d[0], u[0], f, k, term);
  if term = 1 then
    WriteLn('Determinant of A = ', f * intpower(8, k):0:9, ' (k = ', k, ')')
  else
    WriteLn('Incorrect input');

  ReadLn;

end.

Inverse of a matrix (unit inv)

Suppose a square matrix A. Then the matrix A<sup>-1</sup> is the inverse of A if the product A<sup>-1</sup> A is the identity matrix I (all elements zero, except for the diagonal elements which are 1).

The inverse matrix exists if the determinant of A is not zero.

$A^{-1} A = I = \begin{bmatrix} 1 & & 0 \\ & \ddots & \\ 0 & & 1 \end{bmatrix}$

NumLib provides various routines for calculation of the inverse matrix. Which algorithm is used depends on the properties of the matrix A. All routines discussed here have in common that the matrix is stored in the standard NumLib way, i.e. as a 2D array of real numbers with n rows and m columns, or as a 1D array with n m elements.

No special requirements are made for the matrix passed to invgen, except that it is square. The calculation in this routine is based on LU decomposition with partial pivoting. If the matrix is symmetric, it can be processed by invgsy which uses a reduction to tridiagonal form internally. A symmetric positive definite matrix can be passed to invgpd where the inverse is calculated by Cholesky decomposition.

pascal
procedure invgen(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt);   // --> generic matrix
procedure invgsy(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt);   // --> symmetric matrix
procedure invgpd(n, rwidth: ArbInt; var ai: ArbFloat; var term: ArbInt);   // --> symmetric positive definite matrix

Example

Calculate the inverse of the symmetric matrix $A = \left[ \begin{array}{rrrr} 5 & 7 & 6 & 5 \\ 7 & 10 & 8 & 7 \\ 6 & 8 & 10 & 9 \\ 5 & 7 & 9 & 10 \end{array} \right]$ This matrix is symmetric and positive definite (the product b<sup>T</sup> A b with any vector b = [b<sub>1</sub>, b<sub>2</sub>, b<sub>3</sub>]<sup>T</sup> cannot become zero or negative since all elements of this matrix are positive). Therefore, invgpd is best-suited for this task although the other routines can be used as well (uncomment their calls below).

pascal
program inverse_matrix;  

uses
  typ, inv, omv;

const
  n = 4;
  D = 5;

var
  // a is the input matrix to be inverted
  // Note that this is matrix must be symmetric positive definite.
  a: array[1..n, 1..n] of ArbFloat = (
    (5,  7,  6,  5),
    (7, 10,  8,  7),
    (6,  8, 10,  9),
    (5,  7,  9, 10)
  );
  b: array[1..n, 1..n] of Arbfloat;
  c: array[1..n, 1..n] of ArbFloat;
  term: Integer = 0;
  i, j: Integer;

begin
  // write input matrix
  WriteLn('a = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(a[i, j]:10:D);
    WriteLn;
  end;
  WriteLn;

  // Store input matrix because it will be overwritten by the inverse of a
  for i := 1 to n do
    for j := 1 to n do
      b[i, j] := a[i, j];

  // Calculate inverse  -- uncomment the procedure to be used.
  //invgen(n, n, a[1, 1], term);
  //invgsy(n, n, a[1, 1], term);
  invgpd(n, n, a[1, 1], term);

  // write inverse
  WriteLn('a^(-1) = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(a[i, j]:10:D);
    WriteLn;
  end;
  WriteLn;

  // Check validity of result by multiplying inverse with saved input matrix
  omvmmm(a[1, 1], n, n, n,
         b[1, 1], n, n,
         c[1, 1], n);

  // write inverse
  WriteLn('a^(-1) x a = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(c[i, j]:10:D);
    WriteLn;
  end;
end.

Solving systems of linear equations (unit sle)

A system of linear equations (or linear system) is a collection of two or more linear equations involving the same set of variables, x<sub>1</sub> ... x<sub>n</sub>.

$\begin{cases} a_{11} x_1 + a_{12} x_2 + \dots + a_{1n} x_n = b_1 \\ a_{21} x_1 + a_{22} x_2 + \dots + a_{2n} x_n = b_n \\ \ \ \ \vdots \\ a_{m1} x_1 + a_{m2} x_2 + \dots + a_{mn} x_n = b_m \end{cases}$

In order to find a solution which satisfies all equations simultaneously the system of equations is usually converted to a matrix equation A x = b, where

$A = \left [ \begin{array}{ccc} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \dots & a_{mn} \end{array} \right ], \ \mathbf{x} = \left [ \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right ], \ \mathbf{b} = \left [ \begin{array}{c} b_1 \\ b_2 \\ \vdots \\ b_m \end{array} \right ]$

NumLib offers several procedures to solve the matrix equation depending on the properties of the matrix.

Square matrices with standard storage

Here we describe several routines in which the matrix of the equation system is square, i.e. there are as many equations as unknowns. This matrix is assumed to be passed to the procedures in the standard NumLib way as the first element of a 2D or 1D array. The 2D array must be dimensioned to contain at least n rows and n columns, the 1D array must contain at least n<sup>2</sup> elements.

The most general procedure is slegen and applied to any square matrix A. The system is solved using the Gauss elimination procedure with partial pivoting.

If the matrix A is symmetric then the reduction to tridiagonal form generally is more stable. This is done in routine slegsy.

If A, furthermore, is a symmetric positive definite matrix, then the Cholesky method is best-suited. This is applied in routine slegpd.

pascal
procedure slegen(n, rwidth: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt);   // generic matrix
procedure slegsy(n, rwidth: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt);   // symmetric matrix
procedure slegpd(n, rwidth: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt);   // symmetric positive definite matrix

Example

Solve this system of linear equations:

$\begin{cases} 5 x_1 + 7 x_2 + 6 x_3 + 5 x_4 = 57 \\ 7 x_1 + 10 x_2 + 8 x_3 + 7 x_4 = 79 \\ 6 x_1 + 8 x_2 + 10 x_3 + 9 x_4 = 88 \\ 5 x_1 + 7 x_2 + 9 x_3 + 10 x_4 = 86 \end{cases} \qquad \Rightarrow \qquad A= \left[ \begin{array}{rrrr} 5 & 7 & 6 & 5 \\ 7 & 10 & 8 & 7 \\ 6 & 8 & 10 & 9 \\ 5 & 7 & 9 & 10 \end{array} \right] , \; \mathbf{b} = \left[ \begin{array}{r} 57 \\ 79 \\ 88 \\ 86 \end{array} \right]$

This matrix is symmetric positive definite. Therefore, all routines presented here can be applied. Uncomment the requested line in the code below. As a test the result vector x is multiplied to matrix A, and the result must equal the constant vector b.

pascal
program solve_sys_lin_eq;

uses
  typ, sle, omv;

const
  n = 4;
  A: array[1..n, 1..n] of ArbFloat = (
    ( 5,  7,  6,  5),
    ( 7, 10,  8,  7),
    ( 6,  8, 10,  9),
    ( 5,  7,  9, 10)
  );
  b: array[1..n] of ArbFloat = (
    57, 79, 88, 86
  );

var
  x: array[1..n] of ArbFloat;
  b_test: array[1..n] of ArbFloat;
  ca: ArbFloat;
  i, j: Integer;
  term: Integer;

begin
  WriteLn('Solve matrix system A x = b');
  WriteLn;

  // Print A
  WriteLn('Matrix A = ');
  for i:= 1 to n do begin
    for j := 1 to n do
      Write(A[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;

  // Print b
  WriteLn('Vector b = ');
  for i:= 1 to n do
    Write(b[i]:10:0);
  WriteLn;
  WriteLn;

  // Solve - select one of these methods
  slegen(n, n, A[1, 1], b[1], x[1], ca, term);
//  slegsy(n, n, A[1, 1], b[1], x[1], ca, term);
//  slegpd(n, n, A[1, 1], b[1], x[1], ca, term);

  if term = 1 then begin
    WriteLn('Solution vector x = ');
    for i:= 1 to n do
      Write(x[i]:10:0);
    WriteLn;
    WriteLn;

    omvmmv(A[1,1], n, n, n, x[1], b_test[1]);

    WriteLn('Check result: A x = (must be equal to b)');
    for i:= 1 to n do
      Write(b_test[i]:10:0);
    WriteLn;
  end
  else
    WriteLn('Error');

  ReadLn;

end.

Band matrix

A special solution, slegba, is implemented for a band matrix, i.e. a matrix in which all elements are zero outside a band of width l below and of width r above the main diagonal.

pascal
procedure slegba(n, l, r: ArbInt; var a, b, x, ca: ArbFloat; var term:ArbInt);

Example

Solve this system of linear equations:

$\begin{array}{ccccccc} 5 x_1 & - & 4 x_2 & + & x_3 & & & & & & & & & = 0 \\ -4 x_1 & + & 6 x_2 & - & 4 x_3 & + & x_4 & & & & & & & = 0 \\ x_1 & - & 4 x_2 & + & 6 x_3 & - & 4 x_4 & & & & & & & = 0 \\ & & x_2 & - & 4 x_3 & + & 6 x_4 & - & 4 x_5 & & & & & = 1 \\ & & & & x_3 & - & 4 x_4 & + & 6 x_5 & - & 4 x_6 & + & x_7 & = 0 \\ & & & & & & x_4 & - & 4 x_5 & + & 6 x_6 & - & 4 x_7 & = 0 \\ & & & & & & & & x_5 & - & 4 x_6 & + & 5 x_7 & = 0 \\ \end{array} \qquad \Rightarrow \qquad A= \left[ \begin{array}{rrrr} 5 & -4 & 1 & 0 & 0 & 0 & 0 \\ -4 & 6 & -4 & 1 & 0 & 0 & 0 \\ 1 & -4 & 6 & -4 & 1 & 0 & 0 \\ 0 & 1 & -4 & 6 & -4 & 1 & 0 \\ 0 & 0 & 1 & -4 & 6 & -4 & 1 \\ 0 & 0 & 0 & 1 & -4 & 6 & -4 \\ 0 & 0 & 0 & 0 & 1 & -4 & 5 \end{array} \right] , \; \mathbf{b} = \left[ \begin{array}{r} 0 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{array} \right]$

pascal
program solve_band;

uses
  typ, sle, omv;

const
  n = 7;
  L = 2;
  R = 2;
  nelements = n * (L + 1 + R) - (L * (L + 1) + R * (R + 1)) div 2;

  // Band elements arranged in rows, elements outside of the band omitted!
  a: array[1..nelements] of Arbfloat = (
             5, -4, 1,
         -4, 6, -4, 1,
      1, -4, 6, -4, 1,
      1, -4, 6, -4, 1,
      1, -4, 6, -4, 1,
      1, -4, 6, -4,
      1, -4, 5
    );
  b: array[1..n] of ArbFloat = (
    0, 0, 0, 1, 0, 0, 0
  );

var
  x: array[1..n] of ArbFloat;
  a_test: array[1..n, 1..n] of ArbFloat;
  b_test: array[1..n] of ArbFloat;
  ca: ArbFloat;
  i, j, k: Integer;
  term: Integer;

begin
  WriteLn('Solve matrix system A x = b where A is a band matrix');
  WriteLn;

  // Write diagonal elements
  WriteLn('Matrix A:');
  WriteLn('  n = ', n);
  Writeln('  Left band width L = ', L);
  WriteLn('  Right band width R = ', R);
  WriteLn('  Diagonal elements of A = ');
  for k := 1 to nelements do
    Write(a[k]:5:0);
  WriteLn;
  WriteLn;

  // Solve
  slegba(n, l, r, A[1], b[1], x[1], ca, term);

  if term = 1 then begin
    // To test the result we multiply A with the solution vector x and check
    // whether the result is equal to b.
    // Since mutliplication of a matrix stored like a band matrix is quite
    // cumbersome we copy the band matrix into a standard matrix.
    FillChar(a_test, SizeOf(a_test), 0);
    k := 1;
    for i:=1 to L do
      for j := i-L to i+R do
        if (j >= 1) and (j <= n) then begin
          a_test[i,j] := A[k];
          inc(k);
        end;

    for i:= L+1 to n-R do
      for j := i-L to i+R do begin
        a_test[i,j] := A[k];
        inc(k);
      end;
    for i := n-R+1 to n do
      for j := i-L to i+R do
        if j <= n then begin
          a_test[i,j] := A[k];
          inc(k);
        end;

    // Print A
    WriteLn('Matrix A =');
    for i:=1 to n do begin
      for j:=1 to n do
        Write(a_test[i,j]:5:0);
      WriteLn;
    end;
    WriteLn;

    // Print b
    WriteLn('Vector b = ');
    for i:= 1 to n do
      Write(b[i]:10:0);
    WriteLn;
    WriteLn;

    // Print solution
    WriteLn('Solution vector x = ');
    for i:= 1 to n do
      Write(x[i]:10:3);
    WriteLn;
    WriteLn;

    omvmmv(a_test[1,1], n, n, n, x[1], b_test[1]);

    WriteLn('Check result: A x = (must be equal to b)');
    for i:= 1 to n do
      Write(b_test[i]:10:0);
    WriteLn;
  end
  else
    WriteLn('Error');
end.

Symmetric positive definite band matrix

The optimized solution method for a symmetric positive band matrix is called slegpb.

pascal
procedure slegpb(n, w: ArbInt; var a, b, x, ca: ArbFloat; var term: ArbInt);

Example

The matrix of the band matrix example is symmetrical positive definite. Therefore, procedure slegpb can be applied as well:

pascal
program solve_spd;    // "spd" = symmetric positive definite

uses
  typ, sle, omv;

const
  n = 7;
  w = 2;
  nelements = n * (w + 1) - (w * (w + 1)) div 2;

  // Band elements arranged in rows, elements outside of the band as well as right diagonals omitted!
  a: array[1..nelements] of Arbfloat = (
             5,
         -4, 6,
      1, -4, 6,
      1, -4, 6,
      1, -4, 6,
      1, -4, 6,
      1, -4, 5
    );
  b: array[1..n] of ArbFloat = (
    0, 0, 0, 1, 0, 0, 0
  );

var
  x: array[1..n] of ArbFloat;
  a_test: array[1..n, 1..n] of ArbFloat;
  b_test: array[1..n] of ArbFloat;
  ca: ArbFloat;
  i, j, k: Integer;
  term: Integer;

begin
  WriteLn('Solve matrix system A x = b where A is a symmetric positive definite band matrix');
  WriteLn;

  // Write diagonal elements
  WriteLn('Matrix A:');
  WriteLn('  n = ', n);
  Writeln('  (One-sided) band width w = ', w);
  WriteLn('  Diagonal elements of A = ');
  for k := 1 to nelements do
    Write(a[k]:5:0);
  WriteLn;
  WriteLn;

  // Solve
  slegpb(n, w, A[1], b[1], x[1], ca, term);

  if term = 1 then begin
    // To test the result we multiply A with the solution vector x and check
    // whether the result is equal to b.
    // Since mutliplication of a matrix stored like a band matrix is quite
    // cumbersome we copy the band matrix into a standard matrix.
    FillChar(a_test, SizeOf(a_test), 0);
    i := 1;
    k := 1;
    while (k <= nelements) and (i <= n) do begin
      for j := i - w to i do
        if (j >= 1) and (j <= n) then begin
          a_test[i, j] := a[k];
          a_test[j, i] := a[k];
          inc(k);
        end;
      inc(i);
    end;

    // Print A
    WriteLn('Matrix A =');
    for i:=1 to n do begin
      for j:=1 to n do
        Write(a_test[i,j]:5:0);
      WriteLn;
    end;
    WriteLn;

    // Print b
    WriteLn('Vector b = ');
    for i:= 1 to n do
      Write(b[i]:10:0);
    WriteLn;
    WriteLn;

    // Print solution
    WriteLn('Solution vector x = ');
    for i:= 1 to n do
      Write(x[i]:10:3);
    WriteLn;
    WriteLn;

    omvmmv(a_test[1,1], n, n, n, x[1], b_test[1]);

    WriteLn('Check result: A x = (must be equal to b)');
    for i:= 1 to n do
      Write(b_test[i]:10:3);
    WriteLn;
  end
  else
    WriteLn('Error');
end.

Tridiagonal matrix

NumLib offers two specialized procedures for solving a system of linear equations based on a tridiagonal matrix A. Both rely on the memory-optimized storage of tridiagonal matrices as described here. They differ in the internal calculation algorithm and in the behavior with critical matrices.

pascal
procedure sledtr(n: ArbInt; var l, d, u, b, x: ArbFloat; var term: ArbInt);
procedure slegtr(n: ArbInt; var l, d, u, b, x, ca: ArbFloat; var term: ArbInt);

sledtr is numerically stable if matrix A fulfills one of these conditions:

However, this procedure does not provide the parameter ca from which the accuracy of the determined solution can be evaluated. If this is needed the (less stable) procedure slegtr must be used.

Example Solve this tridiagonal system of linear equations for n = 8:

$\begin{cases} 188 x_1 - 100 x_2 = 0 \\ -100 x_1 + 188 x_2 -100 x_3 = 0 \\ \vdots \\ -100 x_{n-2} + 188 x_{n-1} - 100 x_n = 0 \\ -100 x_{n-1} + 188 x_n = 0 \end{cases} \qquad \Rightarrow \qquad A= \left[ \begin{array}{rrrrr} 188 & -100 & & & 0 \\ -100 & 188 & -100 & & \\ & \ddots & \ddots & \ddots & \\ & & -100 & 188 & -100 \\ & & & -100 & 188 \\ \end{array} \right] , \; \mathbf{b} = \left[ \begin{array}{r} 88 \\ -12 \\ \vdots \\ -12 \\ 88 \end{array} \right]$

pascal
program solve_tridiag;

uses
  typ, sle;

const
  n = 8;
  u: array[1..n-1] of ArbFloat = (-100, -100, -100, -100, -100, -100, -100      );
  d: array[1..n] of ArbFloat   = ( 188,  188,  188,  188,  188,  188,  188,  188);
  l: array[2..n] of ArbFloat   = (      -100, -100, -100, -100, -100, -100, -100);
  b: array[1..n] of ArbFloat   = (  88,  -12,  -12,  -12,  -12,  -12,  -12,   88);

var
  A: array[1..n, 1..n] of ArbFloat;
  x: array[1..n] of ArbFloat;
  b_test: array[1..n] of ArbFloat;
  ca: ArbFloat;
  i, j: Integer;
  term: Integer;

begin
  WriteLn('Solve tridiagonal matrix system A x = b');
  WriteLn;

  Write('Superdiagonal of A:':25);
  for i := 1 to n-1 do
    Write(u[i]:10:0);
  WriteLn;

  Write('Main diagonal of A:':25);
  for i:= 1 to n do
    Write(d[i]:10:0);
  WriteLn;

  Write('Subdiagonal of A:':25);
  Write('':10);
  for i:=2 to n do
    Write(l[i]:10:0);
  WriteLn;

  Write('Vector b:':25);
  for i:=1 to n do
    Write(b[i]:10:0);
  WriteLn;

  // Solve
  slegtr(n, l[2], d[1], u[1], b[1], x[1], ca, term);
  { or:
  sledtr(n, l[2], d[1], u[1], b[1], x[1], term);
  }

  if term = 1 then begin
    Write('Solution vector x:':25);
    for i:= 1 to n do
      Write(x[i]:10:0);
    WriteLn;

    // Multiply A with x to test the result
    // NumLib does not have a routine to multiply a tridiagonal matrix with a
    // vector... Let's do it manually.
    b_test[1] := d[1]*x[1] + u[1]*x[2];
    for i := 2 to n-1 do
      b_test[i] := l[i]*x[i-1] + d[i]*x[i] + u[i]*x[i+1];
    b_test[n] := l[n]*x[n-1] + d[n]*x[n];

    Write('Check b = A x:':25);
    for i:= 1 to n do
      Write(b_test[i]:10:0);
    WriteLn;
  end
  else
    WriteLn('Error');
end.

Overdetermined systems (least squares)

Unlike the other routines in unit sle which require a square matrix A, slegls can solve linear systems described by a rectangular matrix which has more rows than colums, or speaking in terms of equations, has more equations than unknowns. Such a system generally is not solvable exactly. But an approximate solution can be found such that the sum of squared residuals of each equation, i.e. the norm $\|\mathbf{b} - A \mathbf{x}\|_2$ is as small as possible. The most prominent application of this technique is fitting of equations to data (regression analysis).

pascal
procedure slegls(var a: ArbFloat; m, n, rwidtha: ArbInt; var b, x: ArbFloat; var term: ArbInt);

The method is based on reduction of the array A to upper triangle shape through Householder transformations.

Example

Find the least-squares solution for the system A x = b of 4 equations and 3 unknowns with

$A= \left( \begin{array}{rrr} 1 & 0 & 1 \\ 1 & 1 & 1 \\ 0 & 1 & 0 \\ 1 & 1 & 0 \end{array} \right), \ \ b= \left( \begin{array}{r} 21 \\ 39 \\ 21 \\ 30 \end{array} \right).$

pascal
program solve_leastsquares;

uses
  typ, sle, omv;

const
  m = 4;
  n = 3;

  A: array[1..m, 1..n] of ArbFloat = (
    (1, 0, 1),
    (1, 1, 1),
    (0, 1, 0),
    (1, 1, 0)
  );
  b: array[1..m] of ArbFloat = (
    21, 39, 21, 30);

var
  x: array[1..n] of ArbFloat;
  term: ArbInt;
  i, j: Integer;
  b_test: array[1..m] of ArbFloat;
  sum: ArbFloat;

begin
  WriteLn('Solve A x = b with the least-squares method');
  WriteLn;

  // Display input data
  WriteLn('A = ');
  for i := 1 to m do
  begin
    for j := 1 to n do
      Write(A[i, j]:10:0);
    WriteLn;
  end;
  WriteLn;
  WriteLn('b = ');
  for i := 1 to m do
    Write(b[i]:10:0);
  WriteLn;

  // Calculate and show solution
  slegls(a[1,1], m, n, n, b[1], x[1], term);

  WriteLn;
  WriteLn('Solution x = ');
  for j:= 1 to n do
    Write(x[j]:10:0);
  WriteLn;

  // Calculate and display residuals
  WriteLn;
  WriteLn('Residuals A x - b = ');
  sum := 0;
  omvmmv(a[1,1], m, n, n, x[1], b_test[1]);
  for i:=1 to m do begin
    Write((b_test[i] - b[i]):10:0);
    sum := sum + sqr(b_test[i] - b[i]);
  end;
  WriteLn;

  // Sum of squared residuals
  WriteLn;
  WriteLn('Sum of squared residuals');
  WriteLn(sum:10:0);

  WriteLn;
  WriteLn('----------------------------------------------------------------------------');
  WriteLn;

  // Modify solution to show that the sum of squared residuals increases';
  WriteLn('Modified solution x'' (to show that it has a larger sum of squared residuals)');
  x[1] := x[1] + 1;
  x[2] := x[2] - 1;
  WriteLn;
  for j:=1 to n do
    Write(x[j]:10:0);
  omvmmv(a[1,1], m, n, n, x[1], b_test[1]);
  sum := 0;
  for i:=1 to m do
    sum := sum + sqr(b_test[i] - b[i]);
  WriteLn;
  WriteLn;
  WriteLn('Sum of squared residuals');
  WriteLn(sum:10:0);
end.

Eigenvalues and eigenvectors (unit eig)

A vector is an eigenvector if it does not change its direction after a linear transformation has been applied to it. In matrix terms: if a non-zero vector x of dimension n is multiplied to a square nxn matrix A and the product is the same vector multiplied by a scalar λ then x is an eigenvector, and λ is called an eigenvalue:

$A \mathbf{x} = \lambda \mathbf{x} \qquad \Rightarrow \qquad A \mathbf{x} - \lambda \mathbf{x} = 0$

NumLib offers a variety of routines for calculation of eigenvectors and eigenvalues. They are optimized for the matrix type. Every matrix type will be supported by two or four procedures of various degree of complexity; they are identified by appended numbers 1 and 2, or 1 to 4:

Matrices with general storage

The routines discussed here assume that the n x n matrix for which the eigenvalues are to be calculated is stored in the conventional way as a 2D (or 1D) array of ArbFloat values.

The NumLib method for determining the eigenvalues and eigenvectors of a generic matrix are called eigge1 and eigge3 ("ge" = "generic"). Optimized calculation schemes exists for generic symmetric matrices and symmetric positive definite matrices; the corresponding routines are called eiggs1, eiggs2, eiggs3, and eiggs4 for generic symmetric matrices ("gs" = "generic symmetric")., and eiggg1, eiggg2, eiggg3, eiggg4 for symmetric positive definite matrices ("gg" = ??)

pascal
// Generic matrix (without any special symmetries)
procedure eigge1(var a: ArbFloat; n, rwidth: ArbInt; var lam: complex; var term: ArbInt);                      // all eigenvalues
procedure eigge3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: complex; rwidthx: ArbInt; var term: ArbInt); // all eigenvalues and eigenvectors

// Generic symmetric matrix
procedure eiggs1(var a: ArbFloat; n, rwidth: ArbInt; var lam: ArbFloat; var term: ArbInt);                     // all eigenvalues
procedure eiggs2(var a: ArbFloat; n, rwidth, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt);             // some eigenvalues (index k1..k2)
procedure eiggs3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: ArbFloat; var term: ArbInt);                 // all eigenvalues and eigenvectors
procedure eiggs4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var lam, x: ArbFloat; var term: ArbInt);         // some eigenvalues and eigenvectors (index k1..k2)

// Symmetric positive definite matrix
procedure eiggg1(var a: ArbFloat; n, rwidth: ArbInt; var lam: ArbFloat; var term: ArbInt);                     // all eigenvalues
procedure eiggg2(var a: ArbFloat; n, rwidth, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt);             // some eigenvalues (index k1..k2)
procedure eiggg3(var a: ArbFloat; n, rwidtha: ArbInt; var lam, x: ArbFloat; var term: ArbInt);                 // all eigenvalues and eigenvectors
procedure eiggg4(var a: ArbFloat; n, rwidtha, k1, k2: ArbInt; var lam, x: ArbFloat; var term: ArbInt);         // some eigenvalues and eigenvectors (index k1..k2)

Additionally, in case of procedure eigge3, eigge4 , eiggg3 or eiggg4:

Additionally, in case of procedures eiggs2, eiggs4, eiggg2 and eiggg4:

Example

Calculate the eigenvalues and eigenvectors of the matrix

$A= \begin{bmatrix} 8 & -1 & -5 \\ -4 & 4 & -2 \\ 18 & -5 & -7 \end{bmatrix}$

pascal
program eig_general_matrix;

uses
  SysUtils, math,
  typ, eig;

const
  n = 3;
  D = 3;

var
  // a is the input matrix
  a: array[1..n, 1..n] of ArbFloat = (
    ( 8, -1, -5),
    (-4,  4, -2),
    (18, -5, -7)
  );
  lambda: array[1..n] of complex;
  x: array[1..n, 1..n] of complex;
  term: Integer = 0;
  i, j: Integer;

  function ComplexToStr(z: complex; Decimals: Integer): String;
  var
    sgn: array[boolean] of string = ('+', '-');
  begin
    Result := Format('%.*f %s %.*fi', [Decimals, z.Re, sgn[z.Im < 0], Decimals, abs(z.Im)]);
  end;

begin
  // write input matrix
  WriteLn('a = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(a[i, j]:10:D);
    WriteLn;
  end;
  WriteLn;

  // Calculate eigenvalues/vectors
  eigge3(a[1,1], n, n, lambda[1], x[1,1], n, term);

  // write eigenvalues
  WriteLn('Eigenvalues: lambda = ');
  for i := 1 to n do
    Write(ComplexToStr(lambda[i], D):25);
  WriteLn;
  WriteLn;

  // Write eigenvectors
  WriteLn('Eigenvectors (as columns): x = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(ComplexToStr(x[i, j], D):25);
    WriteLn;
  end;
end.

Symmetric band matrices

NumLib provides four routines to calculate eigenvalues and eigenvectors of symmetric band matrices: eigbs1...eigbs4 ("bs" = band and symmetric):

pascal
// All eigenvalues
procedure eigbs1(var a: ArbFloat; n, w: ArbInt; var lam: ArbFloat; var term: ArbInt);

// Some eigenvalues (index k1..k2)
procedure eigbs2(var a: ArbFloat; n, w, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt);

// All eigenvalues and all eigenvectors
procedure eigbs3(var a: ArbFloat; n, w: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt; var term: ArbInt);

// Some eigenvalues, some eigenvectors (index k1..k2)
procedure eigbs4(var a: ArbFloat; n, w, k1, k2: ArbInt; var lam, x: ArbFloat; rwidthx: ArbInt; var m2, term: ArbInt);

If the bandwidth w > n/3 then it is better to calculate all eigenvalues and eigenvectors even if not all of them are needed.

Example

Calculate the eigenvalues and eigenvectors of the symmetric 7 x 7 matrix $A= \begin{bmatrix} 5 & -4 & 1 & & & & 0 \\ -4 & 6 & -4 & 1 & & & \\ 1 & -4 & 6 & -4 & 1 & & \\ & 1 & -4 & 6 & -4 & 1 & \\ & & 1 & -4 & 6 & -4 & 1 \\ & & & 1 & -4 & 6 & -4 \\ 0 & & & & 1 & -4 & 5 \end{bmatrix}$ The eigenvalues are $\lambda_k = 16\sin^{4}\frac{k\pi}{16}, \qquad k=1,2,\ldots,7$

pascal
program eig_symmband_matrix;

uses
  SysUtils, math,
  typ, eig;

const
  n = 7;
  w = 2;
  nelements = n * (w + 1) - (w * (w + 1)) div 2;
  D = 3;

  // a contains the elements of the left and main diagonals of the input matrix
  a: array[1..nelements] of ArbFloat = (
     5,
    -4,  6,
     1, -4,  6,
         1, -4,  6,
             1, -4,  6,
                 1, -4, 6,
                     1, -4, 5 );

  function MatrixIndexToArrayIndex(i, j, n, w: Integer): Integer;

    function DiagElementIndex(i: Integer): Integer;
    var
      k: Integer;
    begin
      Result := 1;
      if i = 1 then
        exit;
      // Rows truncated at left
      for k := 2 to w do begin
        Result := Result + k;
        if k = i then
          exit;
      end;
      // full rows and those truncated at right
      for k := w+1 to n do begin
        Result := Result + w + 1;
        if k = i then
          exit;
      end;
      Result := n;
    end;

  var
    d: Integer;
  begin
    if j > i then begin
      Result := MatrixIndexToArrayIndex(j, i, n, w);
      exit;
    end;

    Result := DiagElementIndex(i);
    if (i = j) then
      exit;

    d := i - j;
    if d > w then
      Result := -1
    else begin
      dec(Result, d);
      if (Result < 1) then
        Result := -1
    end;
  end;

var
  lambda: array[1..n] of ArbFloat;
  x: array[1..n, 1..n] of ArbFloat;
  term: Integer = 0;
  i, j, k: Integer;

begin
  // Write diagonal elements
  WriteLn('n = ', n);
  Writeln('(One-sided) band width w = ', w);
  Write('Diagonal elements of A = ', a[1]:0:0);
  for k := 2 to nelements do
    Write(a[k]:3:0);
  WriteLn;
  WriteLn;

  // write reconstructed band input matrix (not needed for calculation)
  WriteLn('Reconstructed A = ');
  k := 1;
  for i := 1 to n do begin
    for j := 1 to n do begin
      k := MatrixIndexToArrayIndex(i, j, n, w);
      if k = -1 then
        Write(0.0:3:0)
      else
        Write(a[k]:3:0);
    end;
    WriteLn;
  end;
  WriteLn;

  // Calculate eigenvalues/vectors
  eigbs3(a[1], n, w, lambda[1], x[1,1], n, term);
  if term <> 1 then begin
    WriteLn('term = ', term, ' --> ERROR');
    halt;
  end;

  // Write expected results of eigenvalues
  WriteLn('Expected eigenvalues:');
  for i := 1 to n do
    Write(16 * intpower(sin(i*pi/16), 4):15:D);
  WriteLn;
  WriteLn;

  // write eigenvalues
  WriteLn('Calculated eigenvalues: lambda = ');
  for i := 1 to n do
    Write(lambda[i]:15:D);
  WriteLn;
  WriteLn;

  // Write eigenvectors
  WriteLn('Eigenvectors (as columns): x = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(x[i, j]:15:D);
    WriteLn;
  end;
end.

Symmetric tridiagonal matrices

This type of matrices is treated best by the routines eigts1..eigts4:

pascal
// All eigenvalues
procedure eigts1(var d, cd: ArbFloat; n: ArbInt; var lam: ArbFloat; var term: ArbInt);

// Some eigenvalues (with indices k1..k2)
procedure eigts2(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam: ArbFloat; var term: ArbInt);

// All eigenvalues and all eigenvectors
procedure eigts3(var d, cd: ArbFloat; n: ArbInt; var lam, x: ArbFloat; rwidth: ArbInt; var term: ArbInt);

// Some eigenvalues and eigenvectors (with indices k1..k2)
procedure eigts4(var d, cd: ArbFloat; n, k1, k2: ArbInt; var lam, x: ArbFloat; rwidth: ArbInt; var m2, term: ArbInt);

Example

Calculate eigenvalues and eigenvectors of the matrix $A= \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 2 & 0 \\ 0 & 2 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{bmatrix}$ The expected eigenvalues are $-\sqrt{2},\ 2-\sqrt{2},\ \sqrt{2},\ 2+\sqrt{2}$

pascal
program eig_symmtridiag_matrix;

uses
  typ, eig;

const
  n = 4;

  // a contains the elements of the main diagonal of the input matrix
  d: array[1..n] of ArbFloat = (1, 1, 1, 1);
  // c contains the elements of the subdiagonal
  c: array[2..n] of ArbFloat = (1, 2, 1);

var
  lambda: array[1..n] of ArbFloat;
  x: array[1..n, 1..n] of ArbFloat;
  term: Integer = 0;
  i, j, k: Integer;

begin
  // Write diagonals elements
  WriteLn('n = ', n);
  Write('Elements of main diagonal = ', d[1]:0:0);
  for k := 2 to n do
    Write(d[k]:3:0);
  WriteLn;
  Write('Elements of subdiagonal   = ', ' ':3, c[2]:0:0);
  for k := 3 to n do
    Write(c[k]:3:0);
  WriteLn;
  WriteLn;

  // write reconstructed band input matrix (not needed for calculation)
  WriteLn('Reconstructed A = ');
  for i := 1 to n do begin
    for j := 1 to n do begin
      if j = i then
        Write(d[i]:3:0)
      else if (j = i-1) then
        Write(c[i]:3:0)
      else if (j = i+1) then
        Write(c[i+1]:3:0)
      else
        Write(0.0:3:0);
    end;
    WriteLn;
  end;
  WriteLn;

  // Calculate eigenvalues/vectors
  eigts3(d[1], c[2], n, lambda[1], x[1,1], n, term);
  if term <> 1 then begin
    WriteLn('term = ', term, ' --> ERROR');
    halt;
  end;

  // Write expected results of eigenvalues
  WriteLn('Expected eigenvalues:');
  Write(-sqrt(2):15:3, 2-sqrt(2):15:3, sqrt(2):15:3, 2+sqrt(2):15:3);
  WriteLn;
  WriteLn;

  // write eigenvalues
  WriteLn('Calculated eigenvalues: lambda = ');
  for i := 1 to n do
    Write(lambda[i]:15:3);
  WriteLn;
  WriteLn;

  // Write eigenvectors
  WriteLn('Eigenvectors (as columns): x = ');
  for i := 1 to n do begin
    for j := 1 to n do
      Write(x[i, j]:15:3);
    WriteLn;
  end;
end.

Finding the roots of a function (unit roo)

The roots are the x values at which a function f(x) is zero.

Roots of a polynomial

A polynomial of degree n

math
z^n + a_1 z^{n-1} + a_2 z^{n-2} + ... + a_{n-1} z + a_n = 0

always has n, not necessarily distinct, complex solutions. The datatype complex has been described in section complex numbers. These roots can be calculated by the function roopol:

pascal
procedure roopol(var a: ArbFloat; n: ArbInt; var z: complex; var k, term: ArbInt);

Example

Calculate the roots of the polynomial z<sup>5</sup> + 3 z<sup>4</sup> + 4 z<sup>3</sup> - 8 z<sup>2</sup>. The expected zero points are:

math
z_1=0,\ z_2=0,\ z_3=1,\ z_4=-2+2i,\ z_5=-2-2i
pascal
program polynomial;

uses
  SysUtils, typ, roo;

const
  n = 5;
  a: array[1..n] of ArbFloat = (3, 4, -8, 0, 0);

var
  z: array[1..n] of complex;
  k: ArbInt;
  term: ArbInt;
  i: Integer;
  c: complex;

  function ComplexToStr(z: complex; Decimals: Integer): string;
  const
    SIGN: array[boolean] of string = ('+', '-');
  begin
    Result := Format('%.*f %s %.*f i', [Decimals, z.re, SIGN[z.im <0], Decimals, abs(z.im)]);
  end;

begin
  // Solve equation
  roopol(a[1], n, z[1], k, term);

  if term = 1 then begin
    // Display results
    WriteLn('Results of procedure roopol:');
    for i:=1 to n do
      WriteLn('  Solution #', i, ': ', ComplexToStr(z[i], 6):20);
    WriteLn;

    // Display expected results
    Writeln('Expected results:');
    c.Init(0, 0);  // z1 = 0
    WriteLn('  Solution #1: ', complexToStr(c, 6):20);
    c.Init(0, 0);  // z2 = 0
    WriteLn('  Solution #2: ', complexToStr(c, 6):20);
    c.Init(1, 0);  // z3 = 1
    WriteLn('  Solution #3: ', complexToStr(c, 6):20);
    c.Init(-2, +2);  // z4 = -2 + 2 i
    WriteLn('  Solution #4: ', complexToStr(c, 6):20);
    c.Init(-2, -2);  // z4 = -2 - 2 i
    WriteLn('  Solution #5: ', complexToStr(c, 6):20);
  end else
    WriteLn('ERROR');
end.

Roots of the quadratic equation

The quadratic equation

${z}^2 + {p} {z} + {q} = 0$

is a special polynomal of degree 2. It always has two, not necessarily different, complex roots. These solutions can be determined by the procedure rooqua.

pascal
procedure rooqua(p, q: ArbFloat; var z1, z2: complex);

Note that rooqua assumes that the coefficient of the quadratic term has been normalized to be 1.

Example

Determine the roots of the equation z<sup>2</sup> + 2 z + 5 = 0.

pascal
program quadratic_equ;

uses
  SysUtils, typ, roo;

var
  z1, z2: complex;

const
  SIGN: array[boolean] of string = ('+', '-');

begin
  rooqua(2, 5, z1, z2);

  WriteLn(Format('1st solution: %g %s %g i', [z1.re, SIGN[z1.im < 0], abs(z1.im)]));
  WriteLn(Format('2nd solution: %g %s %g i', [z2.re, SIGN[z2.im < 0], abs(z2.im)]));
end.

Roots of the binomial equation

The equation

math
z^n = a

is another special polynomial in which all terms except for the highest- and lowest-order terms have zero coefficients. It has n complex solutions which can be calculated by the procedure roobin:

pascal
procedure roobin(n: ArbInt; a: complex; var z: complex; var term: ArbInt);

Example

Calculate the roots of the equation

math
z^4 = -1

The exact solutions are

math
z_1 = \frac{1}{2} \sqrt{2} (1+i) ,\ z_2 = \frac{1}{2} \sqrt{2} (1-i) ,\ z_3= \frac{1}{2} \sqrt{2} (-1+i) ,\ z_4= \frac{1}{2} \sqrt{2} (-1-i)
pascal
program binomial;

uses
  SysUtils, typ, roo;

const
  n = 4;

var
  z: array[1..n] of complex;
  term: ArbInt;
  i: Integer;
  a: complex;
  c: complex;

  function ComplexToStr(z: complex; Decimals: Integer): string;
  const
    SIGN: array[boolean] of string = ('+', '-');
  begin
    Result := Format('%.*g %s %.*g i', [Decimals, z.re, SIGN[z.im < 0], Decimals, abs(z.im)]);
  end;

begin
  // Prepare constant term as a complex value
  a.Init(-1, 0);

  // Solve equation
  roobin(n, a, z[1], term);

  if term = 1 then begin
    // Display results
    WriteLn('Results of procedure roobin:');
    for i:=1 to n do
      WriteLn('  Solution #', i, ': ', ComplexToStr(z[i], 6):20);
    WriteLn;

    // Display expected results
    Writeln('Expected results:');
    c.Init(1, 1);
    c.Scale(0.5*sqrt(2));
    WriteLn('  Solution #1: ', complexToStr(c, 6):20);
    c.Init(1, -1);
    c.Scale(0.5*sqrt(2));
    WriteLn('  Solution #2: ', complexToStr(c, 6):20);
    c.Init(-1, 1);
    c.Scale(0.5*sqrt(2));
    WriteLn('  Solution #3: ', complexToStr(c, 6):20);
    c.Init(-1, -1);
    c.Scale(0.5*sqrt(2));
    WriteLn('  Solution #4: ', complexToStr(c, 6):20);
  end else
    WriteLn('ERROR');
end.

Bisection

In the bisection method two x values a and b are estimated to be around the expected root such that the function values have opposite signs at a and b. The center point of the interval is determined, and the subinterval for which the function has opposite signs at its endpoints is selected for a new iteration. The process ends when a given precision, i.e. interval length, is achieved.

In NumLib, this approach is supported by the procedure roof1r:

pascal
procedure roof1r(f: rfunc1r; a, b, ae, re: ArbFloat; var x: ArbFloat; var term: ArbInt);

Example

The following program determines the square root of 2. This is the x value at which the function f(x) = x^2 - 2 is zero. Since f(1) = 1^2 - 2 = -1 < 0 and f(2) = 2^2 - 2 = 2 > 0 we can assume a and b to be 1 and 2, respectively.

pascal
program bisection_demo;

uses
  typ, roo;

function f(x: ArbFloat): ArbFloat;
begin
  Result := x*x - 2;
end;

var
  x: ArbFloat = 0.0;
  term : ArbInt;

begin
  roof1r(@f, 1.0, 2.0, 1e-9, 0, x, term);
  WriteLn('Bisection result ', x);
  WriteLn('sqrt(2)          ', sqrt(2.0));
end.

Roots of a system of nonlinear equations

Procedure roofnr finds the roots of a system of (nonlinear) equations:

$f_{i}(x_1,x_2,\ldots,x_n)=0, \; i=1,2,\ldots,n$

pascal
procedure roofnr(f: roofnrfunc; n: ArbInt; var x, residu: ArbFloat; ra: ArbFloat; var term: ArbInt);
pascal
  type
    roofnrfunc = procedure(var x, fx: ArbFloat; var deff: boolean);

:* Here, x is the first element of an array of at least n ArbFloat values which provides the x<sub>j</sub> values at which the functions are to be calculated.

:* fx is the first element of another array of n ArbFloat values to receive the result of each function f<sub>i</sub>. Suppose, we want to solve the following system of equations

$\begin{array}{l} f_1(x_1, x_2) = x_1^2 + x_2^2 - 2 = 0 \\ f_2(x_1, x_2) = -(x_1 - 1)^2 + x_2 = 0 \end{array}$

Then the function to be used as f can be written like this:

pascal
    procedure func(var x1, f1x: real; var deff: boolean); 
    var 
      x: array[1..2] of real absolute x1;
      f: array[1..2] of real absolute f1x;
    begin
      f[1] := sqr(x[1]) + sqr(x[2]) - 2;
      f[2] := -sqr(x[1] - 1) + x[2];
    end;

:* The boolean parameter deff provides a way to stop the root finding process if some condition is met. The above example has two solutions, namely the intersection of a circle with a parabola. If only the intersection with x[1] >= 1 is requested then the process can be stopped as soon as a value x[1] < 1 is reached by setting deff to false:

pascal
    procedure func(var x1, f1x: real; var deff: boolean); far;
    var 
      x: array[1..2] of real absolute x1;
      f: array[1..2] of real absolute f1x;
    begin
      deff := x[1] >= 1;
      if deff then begin
        f[1] := sqr(x[1]) + sqr(x[2]) - 2;
        f[2] := -sqr(x[1] - 1) + x[2];
      end;
    end;

$\|f\|_2 = \sqrt{\Sigma_{i=1}^{n} f_i^2}$

Example

Calculate a solutions of the equation system

$\begin{array}{l} f_1(x_1, x_2) = x_1^2 + x_2^2 - 2 = 0 \\ f_2(x_1, x_2) = -(x_1 - 1)^2 + x_2 = 0 \end{array}$

pascal
program nonlinsystem;

uses
  SysUtils, typ, roo;

const
  n = 2;
  ra = 1e-10;

var
  x: array[1..n] of ArbFloat;
  f_check: array[1..n] of ArbFloat;
  residu: ArbFloat;
  i: Integer;
  term: Integer;
  deff: Boolean;

  procedure funcs(var x0, fx: ArbFloat; var deff: boolean);
  var
    xloc: array[1..n] of ArbFloat absolute x0;
    f: array[1..n] of ArbFloat absolute fx;
    x, y, z: ArbFloat;
  begin
    x := xloc[1];
    y := xloc[2];
    f[1] :=  sqr(x) + sqr(y) - 2;
    f[2] := -sqr(x - 1) + y;
  end;

begin
  // Initial guessed values
  x[1] := 0;
  x[2] := 0;

  // Solve equation
  roofnr(@funcs, n, x[1], residu, ra, term);

  WriteLn('term = ', term);
  WriteLn;

  if term in [1, 2] then begin
    // Display results
    WriteLn('Results found by procedure roofnr:');
    for i:=1 to n do
      WriteLn('Solution #' + IntToStr(i)+': ':20, x[i]:0:6);
    WriteLn;

    WriteLn('Norm of residuals: ', residu:0:15);
    WriteLn;

    // Test functions:
    f(x[1], f_check[1], deff);
    WriteLn('Check if the function values really are 0 at the found roots:');
    for i := 1 to n do
      WriteLn('Function #' + IntToStr(i) + ': ':20, f_check[i]:0:15);
  end else
    WriteLn('ERROR');

  ReadLn;
end.

Numerical integration of a function (unit int)

The NumLib function int1fr calculates the integral of a given function between limits a and b with a specified absolute accuracy ae:

pascal
procedure int1fr(f: rfunc1r; a, b, ae: ArbFloat; var integral, err: ArbFloat; var term: ArbInt);

Example

Calculate the integral

$\int_a^b \frac 1 {x^2} \mathrm{d}x$

for several integration limits a and b. Since the function diverges at x = 0 the interval from a to b must not contain this point. The analytic result is $-1/{b} + 1/{a}$

pascal
program integrate;

uses
  SysUtils, typ, int;

function recipx2(x: ArbFloat): ArbFloat;
begin
  Result := 1.0 / sqr(x);
end;

function integral_recipx2(a, b: ArbFloat): Arbfloat;
begin
  if a = 0 then
    a := Infinity
  else if a = Infinity then
    a := 0.0
  else
    a := -1/a;
  if b = 0 then
    b := Infinity
  else if b = Infinity then
    b := 0.0
  else
    b := -1/b;
  Result := b - a;
end;

procedure Execute(a, b: ArbFloat);
var
  err: ArbFloat = 0.0;
  term: ArbInt = 0;
  integral: ArbFloat = 0.0;
begin
  try
    int1fr(@recipx2, a, b, 1e-9, integral, err, term);
  except
    term := 4;
  end;
  Write('  The integral from ' + FloatToStr(a) + ' to ' + FloatToStr(b));
  case term of
    1: WriteLn(' is ', integral:0:9, ', exected: ', integral_recipx2(a, b):0:9);
    2: WriteLn(' is ', integral:0:9, ', error: ', err:0:9, ', exected: ', integral_recipx2(a, b):0:9);
    3: WriteLn(' cannot be calculated: Error in input data');
    4: WriteLn(' cannot be calculated: Divergent, or calculation converges too slowly.');
  end;
end;

begin
  WriteLn('Integral of f(x) = 1/x^2');
  Execute(1.0, 2.0);
  Execute(1.0, 1.0);
  Execute(2.0, 1.0);  
  Execute(1.0, Infinity);
  Execute(-Infinity, -1.0);
  Execute(0.0, Infinity);
  // NOTE: The next line will raise an exception if run in the IDE. This will not happen outside the IDE.
  // Integrate(-1.0, Infinity);
end.

Ordinary differential equations (unit ode)

Solving a single first-order differential equation

The routine odeiv1 solves an initial value problem described by a first-order differential equation of the form $\begin{cases} y' = f(x, y), \qquad x \in [a, b]\\ y(a) = \alpha \end{cases}$ with given a, b, f and initial condition y(a) = α.

pascal
procedure odeiv1(f: rfunc2r; a, ya: ArbFloat; var b, yb: ArbFloat; ae: ArbFloat; var term: ArbInt);

The algorithm is adaptive and is based on an explicit one-step Runge-Kutta method of order five.

The procedure is not suitable for a rigid differential comparison.

The accuracy can not be guaranteed in all cases. In so-called unstable problems, e.g. small variations in y(a) give large variations in y(b), the error may be significantly larger than what the result of the procedure suggests. Call the procedure with different values of ae to investigate this case.

If one wants to solve the initial value problem for a number of points, e.g. from x = 0 to x = 1 with step size 0.1, then it is advisable to "integrate" and thus not to restart at x = 0 with every new step.

Example

Solve the differential equation y" = -10 (y - x<sup>2</sup>) with the initial condition y(0) = 0 for x = 0, 0.5, ..., 5. Compare with the exact solution y(x) = -0.02 exp(-10 x) + x<sup>2</sup> - 0.2 x + 0.02

pascal
program solve_ode;

uses
  typ, ode;

  function f(x, y: ArbFloat): Arbfloat;
  begin
    Result := -10 * (y - sqr(x));
  end;

  function exact(x: real): real; far;
  begin
    Result := -0.02*exp(-10*x) + sqr(x) - 0.2*x + 0.02;
  end;

const
  d = 0.5;         // interval length
  ae = 1e-5;
  n = 10;

var
  a, b: ArbFloat;
  ya, yb: ArbFloat;
  term: ArbInt;
  i: Arbint;

begin
  // Set initial conditions
  a := 0.0;
  b := a + d;
  ya := 0.0;

  WriteLn('x':12, 'y':12, 'exact':12, 'error code':17);
  WriteLn;
  WriteLn(a:12:5, ya:12:5, exact(a):12:5, '-':17);

  for i := 1 to n do begin
    odeiv1(@f, a, ya, b, yb, ae, term);
    WriteLn(b:12:5, yb:12:5, exact(b):12:5, term:17);
    a := b;
    ya := yb;
    b := b + d;
  end;
end.

Solving a system of first-order differential equations

An initial value problem described by system of first-order differential equations can be solved by the procedure odeiv2: $\begin{cases} \mathbf{y}' = \mathbf{f}(x, \mathbf{y}), \qquad x \in [a, b]\\ \mathbf{y}(a) = \mathbf{\alpha} \end{cases}$

where

y: [a, b] --> R<sup>n</sup>

or: y = [y<sub>1</sub>(x), y<sub>2</sub>(x), ..., y<sub>n</sub>(x)]

f: [a, b] --> R<sup>n</sup>

or: f(x, y) = [f<sub>1</sub>(x,y), f<sub>2</sub>(x,y), ..., f<sub>n</sub>(x,y)]

The algorithm is based on an explicit one-step Runge-Kutta method of order five with variable step size.

pascal
procedure odeiv2(f: oderk1n; a: ArbFloat; var ya, b, yb: ArbFloat; n: ArbInt; ae: ArbFloat; var term: ArbInt);

Example

Integrate the following initial value problem between x = 0 and x = 1 with step size 0.1: $\begin{cases} y'_{1}=2xy_{1}+y_{2} \\ y'_{2}=-y_{1}+2xy_{2}\\ y_1(0)=0,\ y_2(0)=1 \end{cases}$

The exact solutions are $\begin{cases} y_1(x)=\exp({x^2}) \sin(x) \\ y_2(x)=\exp({x^2}) \cos(x) \end{cases}$

pascal
program solve_ode_sys;

uses
  typ, ode;

const
  ae = 1e-5;

  procedure f(x: ArbFloat; var y, fxy: ArbFloat);
  var
    _y: array[1..2] of ArbFloat absolute y;
    _fxy: array[1..2] of ArbFloat absolute fxy;
  begin
    _fxy[1] := 2*x*_y[1] + _y[2];
    _fxy[2] := -_y[1] + 2*x*_y[2]
  end;

  function exact1(x: ArbFloat): ArbFloat;
  begin
    Result := exp(x*x) * sin(x);
  end;

  function exact2(x: ArbFloat): ArbFloat;
  begin
    Result := exp(x*x) * cos(x);
  end;

var
  a, b, d: ArbFloat;
  ya: array[1..2] of ArbFloat;
  yb: array[1..2] of ArbFloat;
  term: ArbInt;
  i: Arbint;
  n: ArbInt;

begin
  // Set initial conditions
  a := 0.0;
  b := 0.1;
  d := b - a;
  ya[1] := 0.0;
  ya[2] := 1.0;
  n := 10;

  WriteLn('x':12, 'y[1]':12, 'y[2]':12, 'exact[1]':12, 'exact[2]':12, 'error code':17);
  WriteLn;
  WriteLn(a:12:5, ya[1]:12:5, ya[2]:12:5, exact1(a):12:5, exact2(a):12:5, '-':17);

  for i := 1 to n do begin
    odeiv2(@f, a, ya[1], b, yb[1], 2, ae, term);
    WriteLn(b:12:5, yb[1]:12:5, yb[2]:12:5, exact1(b):12:5, exact2(b):12:5, term:17);
    a := b;
    ya[1] := yb[1];
    ya[2] := yb[2];
    b := b + d;
  end;
end.

Interpolation and fitting (unit ipf)

Unit ipf contains routines for

Fitting means: determine an approximating smooth function such that the deviations from the data points are at minimum.
Interpolation means: determine a function that runs exactly through the data points.

Use of polynomials or splines is recommended unless the data are known to belong to a known function in which there are still some unknown parameters, for example, the data are measurements of the function ''f(x) = a + b e<sup>-c x</sup>.

In the most common case, the data contain measurement errors and appear to follow the shape of a "simple" function. Therefore, it is recommended to first try to fit with a low-degree polynomial, e.g. not higher than degree 5. If this is not successful or if the shape of the function is more complicated try to fit with a spline.

If the data contain no (or very little) measurement noise, it is almost always wise to interpolate the data with a spline.

Fitting with a polynomial

For this purpose, the routine ipfpol can be applied. It uses given data points (x<sub>i</sub>, y<sub>i</sub>), i = 1, ..., m, and n given coefficients b<sub>0</sub>, ..., b<sub>n</sub> to calculate a polynomial $\operatorname{p}(x) = b_0 + b_1 x + \dots + b_n x^n$

In the sense of the least squares fitting, the coefficients b<sub>0</sub> ... b<sub>n</sub> are adjusted to minimize the expression

$\sum_{i=1}^m {(\operatorname{p}(x) - y_i)^2}$

pascal
procedure ipfpol(m, n: ArbInt; var x, y, b: ArbFloat; var term: ArbInt);

After the coefficients b<sub>0</sub>, ..., b<sub>n</sub> once have been determined the value of p(x) can be calculated for any x using the procedure spepol from the unit spe.

Note: With this routine, of course, an interpolating polynomial can also be determined by equaling n to m - 1.

Example

Fit a polynomial of degree 2 to these data points: $`\begin{array}{ccc} \hline i & x_i & y_i \ \hline

pascal
1 & 0.00 & 1.26 \\
2 & 0.08 & 1.37 \\
3 & 0.22 & 1.72 \\
4 & 0.33 & 2.08 \\
5 & 0.46 & 2.31 \\
6 & 0.52 & 2.64 \\
7 & 0.67 & 3.12 \\
8 & 0.81 & 3.48

\end{array}`$

pascal
program polynomial_fit;

uses
  SysUtils, StrUtils, typ, ipf, spe;

const
  m = 8;
  n = 2;
  x: array[1..m] of Arbfloat = (
    0.00, 0.08, 0.22, 0.33, 0.46, 0.52, 0.67, 0.81
  );
  y: array[1..m] of Arbfloat = (
    1.26, 1.37, 1.72, 2.08, 2.31, 2.64, 3.12, 3.48
  );

var
  b: array[0..n] of ArbFloat;
  i: Integer;
  term: ArbInt;
  xint, yint: ArbFloat;
begin
  WriteLn('Fitting a polynomial of degree 2: p(x) = b[0] + b[1] x + b[2] x^2');
  WriteLn;

  WriteLn('Data points');
  WriteLn('i':10, 'xi':10, 'yi':10);
  for i:=1 to m do
    WriteLn(i:10, x[i]:10:2, y[i]:10:2);
  WriteLn;

  // Execute fit
  ipfpol(m, n, x[1], y[1], b[0], term);

  // Display fitted polynomial coefficients
  WriteLn('Fitted coefficients b = ');
  WriteLn(b[0]:10:6, b[1]:10:6, b[2]:10:6);
  WriteLn;

  // Interpolate and display fitted polynomial for some x values
  WriteLn('Some interpolated data');
  WriteLn('':10, 'x':10, 'y':10);
  for i:= 1 to 5 do begin
    xint := 0.2*i;
    yint := spepol(xint, b[0], n);
    WriteLn('':10, xint:10:2, yint:10:2);
  end;
end.

Interpolation with a natural cubic spline

The routine ipfisn can be used for this purpose. It calculates the parameters of the spline. Once these parameters are known, the value of the spline can be calculated at any point using the ipfspn procedure.

Assume that (x<sub>i</sub>, y<sub>i</sub>*),i= 0, ...,*n' are the given data points with x<sub>0</sub> < x<sub>1</sub> < ... < x<sub>n</sub>. Then assume a function g(x) which is twice continuously differentiable (i.e., the first and second derivatives exist and are continuous) and which interpolates the data (i.e., g(x<sub>i</sub> = y<sub>i</sub>). Imagine g(x) as the shape of a thin flexible bar. Then the bending energy (in first approximation) is proportional to

$\int_{-\infty}^{\infty} g''(x)^2 dx$

The interpolating natural spline s can be understood as the function g which minimizes this integral and therefore the bending energy.

The spline s has the following properties:

pascal
procedure ipfisn(n: ArbInt; var x, y, d2s: ArbFloat; var term: ArbInt);
pascal
function  ipfspn(n: ArbInt; var x, y, d2s: ArbFloat; t: ArbFloat; var term: ArbInt): ArbFloat;

Example

Interpolate a natural cubic spline through the following data points: $`\begin{array}{ccc} \hline i & x_i & y_i \ \hline

pascal
0 & 0.00 & 0.980 \\
1 & 0.09 & 0.927 \\
2 & 0.22 & 0.732 \\
3 & 0.34 & 0.542 \\
4 & 0.47 & 0.385 \\
5 & 0.58 & 0.292 \\
6 & 0.65 & 0.248 \\
7 & 0.80 & 0.179 \\
8 & 0.93 & 0.139

\end{array}`$

pascal
program spline_interpolation;

uses
  typ, ipf;

const
  n = 8;
  x: array[0..n] of Arbfloat = (
    0.00, 0.09, 0.22, 0.34, 0.47, 0.58, 0.65, 0.80, 0.93
  );
  y: array[0..n] of Arbfloat = (
    0.990, 0.927, 0.734, 0.542, 0.388, 0.292, 0.248, 0.179, 0.139
  );

var
  s: array[0..n] of ArbFloat;
  d2s: array[1..n-1] of ArbFloat;
  i: Integer;
  term: ArbInt;
  xint, yint: ArbFloat;
begin
  WriteLn('Interpolation with a natural cubic spline');
  WriteLn;

  WriteLn('Data points');
  WriteLn('i':10, 'xi':10, 'yi':10);
  for i:=0 to n do
    WriteLn(i:10, x[i]:10:2, y[i]:10:3);
  WriteLn;

  // Interpolation
  ipfisn(n, x[0], y[0], d2s[1], term);

  // Display 2nd derivatives of splines at xi
  WriteLn('s"(xi)');
  Write(0.0:10:6);      // 2nd derivative of 1st point is 0
  for i:=1 to n-1 do
    Write(d2s[i]:10:6);
  WriteLn(0.0:10:6);    // 2nd derivative of last point is 0
  WriteLn;

  // Interpolate for some x values
  WriteLn('Some interpolated data');
  WriteLn('':10, 'x':10, 'y':10);
  for i:= 1 to 5 do begin
    xint := 0.2 * i;
    yint := ipfspn(n, x[0], y[0], d2s[1], xint, term);
    WriteLn('':10, xint:10:2, yint:10:2);
  end;

end.

Fitting with a natural cubic spline

The routine ipffsn can be used for this purpose. This routine calculates the spline-fit parameters. Once these parameters are known the value of the spline can be calculated at any point using the procedure ipfspn.

Following the model with the bending energy in the previous chapter the calculated fit minimizes

$\int_{-\infty}^{\infty} g''(x)^2 dx + \lambda \sum_{i=0}^n (g(x_i) - y_i)^2$

where the parameter λ must be determined. It is a measure of the compromise between the "smoothness" of the fit - which is measured by the integral (the bending energy) - and the deviation of the fit to the given data - which is measured by the sum of squares. The solution is again a natural cubic spline.

Warning: The procedure ipffsn is buggy and crashes. Therefore, it is not documented at the moment.

Special functions (unit spe)

Evaluation of a polynomial

Horner's scheme provides an efficient method for evaluating a polynomial at a specific x value:

$\operatorname{p}(x) = a_0 + a_1 x + a_2 x^2 + ... + a_n x^n \ \Rightarrow \ \operatorname{p}(x) = a_0 + x (a_1 + x (a_2 + \dots + x (a_{n-1} + x a_n)))$

This method is applied in NumLib's procedure spepol:

pascal
function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;

Example

Evaluate the polynomial

$\operatorname{p}(x) = 2 x^3 - 6 x^2 + 2 x - 1$

at x = 3. Note that the terms are not in the correct order for spepol and rearrange them to

$\operatorname{p}(x) = -1 + 2 x - 6 x^2 + 2 x^3$

pascal
program polynomial;

uses
  typ, spe;

const
  n = 3;
  a: array[0..n] of ArbFloat = (-1, 2, -6, 2);

var
  x, y: ArbFloat;
begin
  x := 3.0;
  y := spepol(x, a[0], n);
  WriteLn(y:0:6);
end.

Error function

The error function erf(x) and its complement, the complementary error function erfc(x), are the lower and upper integrals of the Gaussian function, normalized to unity:

$\operatorname{erf}(x) = \frac {2} {\sqrt{\pi}} \int_0^x e^{-{t}^2} dt$

$\operatorname{erfc}(x) = \frac {2} {\sqrt{\pi}} \int_x^\infty e^{-{t}^2} dt = {1} - \operatorname{erf}(x)$

Both functions can be calculated by the NumLib functions speerf and speefc, respectively

pascal
function speerf(x: ArbFloat): ArbFloat;   // --> erf(x)
function speefc(x: ArbFloat): ArbFloat;   // --> erfc(x)

Example

A table of some values of the error function and the complementary error function is calculated in the following demo project:

pascal
program erf_Table;

uses
  SysUtils, StrUtils, typ, spe;

const
  Wx = 7;
  Wy = 20;
  D = 6;

var
  i: Integer;
  x: ArbFloat;
  fs: TFormatSettings;
begin
  fs := DefaultFormatSettings;
  fs.DecimalSeparator := '.';
  WriteLn('x':Wx, 'erf(x)':Wy+1, 'erfc(x)':Wy+1);
  WriteLn;
  for i:= 0 to 20 do begin
    x := 0.1*i;
    WriteLn(Format('%*.1f %*.*f %*.*f', [
      Wx, x,
      Wy, D, speerf(x),
      Wy, D, speefc(x)
    ], fs));
  end;
end.

Normal distribution

The (cumulative) normal distribution is the integral of the Gaussian function with zero mean and unit standard deviation, i.e.

$\operatorname{N}(x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{x} \operatorname{exp}(\frac{t^2}{2}) dt = \frac{1}{2}[1 + \operatorname{erf} (\frac{x}{\sqrt{2}})]$

Suppose a set of random floating point numbers which are distributed according to this formula. Then N(x) measures the probability that a number from this set is smaller than x.

In NumLib, the normal distribution is calculated by the function normaldist:

pascal
function normaldist(x: ArbFloat): ArbFloat;

The inverse normal distribution of y returns the x value for which the normal distribution N(x) has the value y.

pascal
function invnormaldist(y: ArbFloat): ArbFloat;

Note: The functions normaldist and invnormaldist are not available in fpc 3.0.x and older.

Example

Calculate the probabilities for x = -10, -8, ..., 12 for a normal distribution with mean μ = 2 and standard deviation σ = 5. In order to reduce this distribution to the standard normal distribution implemented by NumLib, we introduce a transformation x<sub>r</sub> = (x - μ)/σ.

pascal
program normdist;

uses
  typ, spe;

const
  mu = 2.0;
  sigma = 5.0;
  x: array[0..12] of ArbFloat = (-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14);

var
  xr, y: ArbFloat;
  i: Integer;
begin
  WriteLn('Normal distribution with mean ', mu:0:1, ' and standard deviation ', sigma:0:1);
  WriteLn;
  WriteLn('x':20, 'y':20);

  for i:= 0 to High(x) do begin
    xr := (x[i] - mu) / sigma;
    y := normaldist(xr);
    WriteLn(x[i]:20:1, y:20:6);
  end;
end.

Gamma function

The gamma function is needed by many probability functions. It is defined by the integral

$\Gamma({x}) = \int_0^{\infty}t^{x-1} e^{-t} dt$

NumLib provides two functions for its calculation:

pascal
function spegam(x: ArbFloat): ArbFloat;
function spelga(x: ArbFloat): ArbFloat;

The first one, spegam, calculates the function directly. But since the gamma function grows rapidly for even not-too large arguments this calculation very easily overflows.

The second function, spelga calculates the natural logarithm of the gamma function which is more suitable to combinatorial calculations where multiplying and dividing the large gamma values can be avoided by adding or subtracting their logarithms.

Example

The following demo project prints a table of some values of the gamma function and of its logarithm. Note that spegam(x) overflows above about x = 170 for data type extended.

pascal
program Gamma_Table;

uses
  SysUtils, StrUtils, typ, spe;

const
  VALUES: array[0..2] of ArbFloat = (1, 2, 5);
  Wx = 7;
  Wy = 30;
  Wln = 20;

var
  i: Integer;
  x: ArbFloat;
  magnitude: ArbFloat;
begin
  WriteLn('x':Wx, 'Gamma(x)':Wy, 'ln(Gamma(x))':Wln);
  WriteLn;
  magnitude := 1E-3;
  while magnitude <= 1000 do begin
    for i := 0 to High(VALUES) do begin
      x := VALUES[i] * magnitude;
      Write(PadLeft(FloatToStr(x), Wx));
      if x <= 170 then  // Extended overflow above 170
        Write(FormatFloat('0.000', spegam(x)):Wy)
      else
        Write('overflow':Wy);
      WriteLn(spelga(x):Wln:3);
      if abs(x-1000) < 1E-6 then
        break;
    end;
    magnitude := magnitude * 10;
  end;
end.

Incomplete gamma function

The lower and upper incomplete gamma functions are defined by

$\operatorname{P}({s},{x}) = \frac{1}{\Gamma({s})} \int_0^{x}t^{s-1} e^{-t} dt$

$\operatorname{Q}({s},{x}) = \frac{1}{\Gamma({s})} \int_{x}^{\infty}t^{s-1} e^{-t} dt = 1 - \operatorname{P}({s}, {x})$

pascal
function gammap(s, x: ArbFloat): ArbFloat;
function gammaq(s, x: ArbFloat): ArbFloat;

Note: The functions gammap and gammaq are not available in fpc v3.0.x and older.

Example

The following project prints out a table of gammap and gammaq values:

pascal
program IncompleteGamma_Table;

uses
  SysUtils, StrUtils, typ, spe;

const
  s = 1.25;

  Wx = 7;
  Wy = 20;
  D = 6;

var
  i: Integer;
  x: ArbFloat;
  fs: TFormatSettings;
begin
  fs := DefaultFormatSettings;
  fs.DecimalSeparator := '.';
  WriteLn('x':Wx, 'GammaP('+FloatToStr(s,fs)+', x)':Wy+1, 'GammaQ('+FloatToStr(s,fs)+', x)':Wy+1);
  WriteLn;
  for i := 0 to 20 do begin
    x := 0.5*i;
    WriteLn(Format('%*.1f %*.*f %*.*f', [
      Wx, x,
      Wy, D, gammap(s, x),
      Wy, D, gammaq(s, x)
    ], fs));
  end;
end.

Beta function

The Beta function is another basis of important functions used in statistics. It is defined as

$\operatorname{B}(a, b) = \frac{\Gamma(a+b)} = \int_0^1{t^{x-1} (1-t)^{y-1} dt}$

The NumLib routine for its calculation is

pascal
function beta(a, b: ArbFloat): ArbFloat;

Note: The function beta is not available in fpc v3.0.x or older.

Incomplete beta function

The (regularized) incomplete beta function is defined by

$\operatorname{I}_x(a,b) = \frac {1}{\operatorname{B}(ab)} \int_0^x{t^{x-1} (1-t)^{y-1} dt}$

pascal
function betai(a, b, x: ArbFloat): ArbFloat;

The inverse incomplete beta function determines the argument x for which the passed parameter is equal to I<sub>x</sub>(a, b):

pascal
function invbetai(a, b, y: ArbFloat): ArbFloat;

Note: The functions betai and invbetai are not available in fpc v3.0.x and older.

Example

Calculate the incomplete beta function for some values of x, a, b, and apply the result to the inverse incomplete beta function to demonstrate the validity of the calculations.

pascal
program incomplete_beta_func;

uses
  SysUtils, typ, spe;

  procedure TestBetai(a, b, x: ArbFloat);
  var
    xc, y: ArbFloat;
  begin
    y := betai(a, b, x);
    xc := invbetai(a, b, y, macheps);
    WriteLn('  x = ', x:0:6, '  --->  y = ', y:0:9, '  --->  x = ', xc:0:6);
  end;

begin
  WriteLn('Incomplete beta function and its inverse');
  WriteLn;
  WriteLn('        y = betai(1, 3, x)      x = invbetai(1, 3, y)');
  TestBetaI(1, 3, 0.001);
  TestBetaI(1, 3, 0.01);
  TestBetaI(1, 3, 0.1);
  TestBetaI(1, 3, 0.2);
  TestBetaI(1, 3, 0.9);
  TestBetaI(1, 3, 0.99);
  TestBetaI(1, 3, 0.999);
  WriteLn;

  WriteLn('        y = betai(1, 1, x)      x = invbetai(1, 1, y)');
  TestBetaI(1, 1, 0.001);
  TestBetaI(1, 1, 0.01);
  TestBetaI(1, 1, 0.1);
  TestBetaI(1, 1, 0.2);
  TestBetaI(1, 1, 0.9);
  TestBetaI(1, 1, 0.99);
  TestBetaI(1, 1, 0.999);
  WriteLn;

  WriteLn('        y = betai(3, 1, x)      x = invbetai(3, 1, y)');
  TestBetaI(3, 1, 0.001);
  TestBetaI(3, 1, 0.01);
  TestBetaI(3, 1, 0.1);
  TestBetaI(3, 1, 0.2);
  TestBetaI(3, 1, 0.9);
  TestBetaI(3, 1, 0.99);
  TestBetaI(3, 1, 0.999);
  WriteLn;
end.

Bessel functions

The Bessel functions are solutions of the Bessel differential equation:

${x}^2 y'' + x y' + ({x}^2 - \alpha^2) {y} = 0$

The Bessel functions of the first kind, or cylindrical harmonics, J<sub>α</sub>(x), are the solutions which are not singular at the origin, in contrast to the Bessel functions of the second kind, Y<sub>α</sub>(x). Y<sub>α</sub>(x) is only defined for x > 0.

Example

The solutions for a purely imaginary argument are called modified Bessel functions of the first kind, I<sub>α</sub>(x), and modified Bessel function of the second kind, K<sub>α</sub>(x). K<sub>α</sub>(x) is only defined for x > 0.

NumLib implements only the solutions for the parameters α = 0 and α = 1:

pascal
function spebj0(x: ArbFloat): ArbFloat;     // Bessel function of the first kind, J0 (alpha = 0)
function spebj1(x: ArbFloat): ArbFloat;     // Bessel function of the first kind, J1 (alpha = 1)

function speby0(x: ArbFloat): ArbFloat;     // Bessel function of the second kind, Y0 (alpha = 0)
function speby1(x: ArbFloat): ArbFloat;     // Bessel function of the second kind, Y1 (alpha = 1)

function spebi0(x: ArbFloat): ArbFloat;     // modified Bessel function of the first kind, I0 (alpha = 0)
function spebi1(x: ArbFloat): ArbFloat;     // modified Bessel function of the first kind, I1 (alpha = 1)

function spebk0(x: ArbFloat): ArbFloat;     // modified Bessel function of the second kind, K0 (alpha = 0)
function spebk1(x: ArbFloat): ArbFloat;     // modified Bessel function of the second kind, K1 (alpha = 1)

Other functions

The unit spe contains some other functions which will not be documented here because they are already available in th FPC's standard unit math.

<thead>

in <em>math</em></th>

</thead> <tbody>

</tbody>

From the Free Pascal wiki, page "NumLib Documentation", as it stood on 2024-02-14 after a change by Wp, out of the wiki's dump of 2024-05-26. Copyright stays with the author who wrote it, according to the old wiki's Licensing page. If you have concerns about this material, message the admin.

Revision 2, last changed by WikiArchive .

Talk about this page

Nobody has opened a bench for this page yet. Sign in to open one.