## Scientific programming with Delphi/Kylix

 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:

• Unit creation
• Components: Memo, OpenDialog, SaveDialog, Panel and Button
• Local functions/procedures
• Text File I/O
• Try...finally statement

#### 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])+', ';
end;
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
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;
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.

(in preparation)

(in preparation)

### Enhancement

• Use of  IPolynom and IPolySystem (TPolynom and TPolySystem) from my PasMatLib

### Other projects

 Scientific Programming