Scientific programming with Delphi/Kylix

Up Scientific Programming

Orthogonal Polynomials project

This program calculates the set of polynomials orthogonal on a discrete set of points x1, . . , xM . This is a Pascal implementation of the Orthogonal basis functions tutorial.

Subject covered:

Algorithm

The initial normalized polynomial is 
 

The polynomials are built recursively: 

   (non-normalized)
     (normalized)

In other words, we start with monomials (polynomials with exactly one nonzero summand)

, n = 0, 1, . . . , M-1 

and transform them to the orthogonal and normalized set of polynomials. 

Design

Now we separate numerical code and the input/output interface. 

Start Delphi and save Unit1 as fPolyOrt and Project1 as PolyOrt in SciProg/PolyOrt directory. Go to File/New/Unit (or  File/New/Other/Unit, depending on the version). Save new Unit1 as uPolyOrt. Add uPolyOrt to the uses clause of the fPolyOrt unit. Try to compile and check for errors.

Visual Design

Drop a Panel on the form (Align = alLeft), the a Memo (Align = alClient), an OpenDialog and a SaveDialog. 

Select the panel and drop 4 buttons: ButtonLoad, ButtonMono, ButtonGenerate and ButtonSave.

Numerical code (uPolyOrt.pas)

 Delphi has two native types of arrays:

  1. Static arrays with arbitrary boundaries, for example, S: array[1 . . 5, -2 . . 2] of string
  2. Dynamic zero-based arrays, for example, B: array of array of integer - 2D integer array

Delphi allows creating dynamic array objects with arbitrary boundaries (see, for example, my PasMatLib library). 

In this tutorial we use the second type of arrays. Zero-based arrays are naturally easier to implement for our problem because the set of polynomials starts with the zero degree.

In the interface part of the unit put the following type, variable and procedure declaration:

unit uPolyOrt;

interface

type
  TFloat = double; // can be changed in one place to another type, e.g. extended

  TInt = integer;

  TPolynomials = array of array of TFloat;  // dynamic 2D array (matrix)

  TCoordinates = array of TFloat;           // dynamic 1D array (vector)

var
// Global variable declarations

  Polynomials: TPolynomials;

  Coord: TCoordinates;

// procedure declarations

  { Sets length for vector X[0..n-1]
    and makes a triangular 2d array
    X
    X X
    X X X
    X X X X
    . . . . . }
  procedure InitArrays(n: TInt; var  X: TCoordinates; var P: TPolynomials);

  { Pn(x) = x**n }
  procedure SetMonomials(P: TPolynomials);

  procedure Orthonormalize(P: TPolynomials; X: TCoordinates);

In the type declaration section we define aliases for several standard types. There are two major reasons for this. It is easy to re-declare the precision (ie. single, double, extended) and make type names either shorter or more descriptive.

Type the following implementation code. 

implementation

// procedure implementations

procedure InitArrays(n: TInt; var X: TCoordinates; var P: TPolynomials);
var
  i: TInt;
begin
  SetLength(X,n);
  SetLength(P,n);
  for i := Low(P) to High(P) do   // triangular matrix
    SetLength(P[i], i+1);
end;

procedure SetMonomials(P: TPolynomials);
var
  i1, i2: TInt;
begin
  for i1 := Low(P) to High(P) do
  begin
    for i2 := Low(P[i1]) to High(P[i1]) do
      P[i1,i2] := 0.0;
    P[i1,i1] := 1.0;
  end;
end;

procedure Orthonormalize(P: TPolynomials; X: TCoordinates);
var
  i, j, k, M: TInt;
  f: TFloat;

  // internal functions local to the procedures Orthonormalize
  
  // Eval(i,x) = Pi(x)
  function Eval(i: TInt; x: TFloat): TFloat;
  var
    j: TInt;
    s: TFloat;
  begin
    s := 0.0;
    for j := High(P[i]) downto Low(P[i]) do
    begin
      s := P[i,j] + x*s;  // x - global to Eval
    end;
    result := s;
  end;

  // Dot product of two polynomials: DotProd(i, j) = <Pi(x)|Pj(x)>
  function DotProd(i, j: TInt): TFloat;
  var
    k: TInt;
    s: TFloat;
  begin
    s := 0.0;
    for k := Low(X) to High(X) do
    begin
      s := s + Eval(i, X[k])*Eval(j, X[k]) // X - global to DotProd
    end;
    result := s;
  end;

begin
  M := Length(X);
  f := M;

  P[0,0] := 1.0/sqrt(f);   // P0

  for i := 1 to M-1 do     // P1, P2, ...
  begin
    for j := 0 to i-1 do   // subtract all previous from Pi
    begin
      f := DotProd(i, j);                 //  <Pi(x)|Pj(x)>
      for k := Low(P[j]) to High(P[j]) do // subtract Pj
        P[i,k] := P[i,k] - f*P[j,k];
    f := DotProd(i, i);                   //  <Pi(x)|Pi(x)>
    for k := Low(P[i]) to High(P[i]) do   // normalize
      P[i,k] := P[i,k]/sqrt(f);
    end;
  end;
end;

end.

Delphi dynamics arrays allow the creation of triangular matrices. However the algorithm works fine with square matrices as well. For the use of square matrices the code  

  SetLength(P,n);
  for  i := Low(P) to High(P) do SetLength(P[i], i+1);

in InitArrays should be replaced with

  SetLength(P,n,n);

You can open and view Polynomial.dpr (program unit) and fPolynomial.pas (form unit). Delphi also generates several files with IDE settings and resources (see "Delphi generated files" in Delphi Help).

User interface and Input/Output (fPolyOrt.dpr)

 Add the following declaration to the private section of TForm1:

    procedure DisplayP(P: TPolynomials);

Press Ctrl-SHIFT-c and Delphi creates the empty procedure in the implementation section. Insert the following code:

procedure TForm1.DisplayP(P: TPolynomials);
var
  i1, i2: TInt;
  s: string;
begin
  for i1 := Low(P) to High(P) do
  begin
    s := '';
    for i2 := Low(P[i1]) to High(P[i1]) do
      s := s + FloatToStr(P[i1,i2])+', ';
    Memo1.Lines.Add(s);
  end;
  Memo1.Lines.Add('');
end;

After double-clicking on each button, add the appropriate code into the generated event handler templates.

procedure TForm1.ButtonLoadClick(Sender: TObject);
var
  f: TextFile;
  k, n: TInt;
  t: TFloat;
  s: string;
begin
  if OpenDialog1.Execute then
  try
    AssignFile(f,OpenDialog1.FileName);
    Reset(f); // Open existing file
    ReadLn(f,n);
    InitArrays(n, Coord, Polynomials);
    Memo1.Clear;
    s := ''; // empty string;
    for k := Low(Coord) to High(Coord) do
    begin
      ReadLn(f,t);   // one coordinate per line
      Coord[k] := t;
      s := s + FloatToStr(t)+', ';
    end;
  finally
    CloseFile(f) // close file regardless of errors
  end;
  Memo1.Lines.Add(s);
  Memo1.Lines.Add('');
end;

procedure TForm1.ButtonMonoClick(Sender: TObject);
begin
  SetMonomials(Polynomials);
  DisplayP(Polynomials);
end;

procedure TForm1.ButtonGenerateClick(Sender: TObject);
begin
  Orthonormalize(Polynomials, Coord);
  DisplayP(Polynomials);
end;

{ Formatted output of comma-separated coefficients }
procedure TForm1.ButtonSaveClick(Sender: TObject);
var
  i1, i2: TInt;
  s: string;
  f: TextFile;
begin
  if SaveDialog1.Execute then
  try
    AssignFile(f,SaveDialog1.FileName);
    Rewrite(f); // Open new or reset existing file
    s := ''; // empty string;
    for i1 := Low(Polynomials) to High(Polynomials) do
    begin
      s := '';
      for i2 := Low(Polynomials[i1]) to High(Polynomials[i1])-1 do
        Write(f,Polynomials[i1,i2]:15:12,',');
      WriteLn(f,Polynomials[i1,i1]:15:12);  // the lasr coeff and EndLine
    end;
  finally
    CloseFile(f) // close file regardless of errors
  end;
end;

Data file preparation

To reproduce the result from Orthogonal basis functions tutorial: Go to   File/New/Other/Text. Type :

 4
 -0.5
 -0.16666666666667
 0.16666666666667
 0.5

Save the file as coord.txt

Test run

Compile and run the program. Press the "Load" button and select coord.txt. Press the "Set Monomials" button and then the "Generate" button. Your results should coincide with the example from the Orthogonal basis functions tutorial within the round-off errors. You can save calculated coefficients in the file of your choice by pressing the "Save" button. 

 The project code for Delphi 7 can be downloaded here.

Kylix-specific notes

(in preparation)

Lazarus-specific notes

(in preparation)

Enhancement 

Other projects

Up Scientific Programming

Home | Resumé |  Shokhirev.com |  Computing |  Links |  Publications

©Nikolai V. Shokhirev, 2001-2005