|
Scientific Programming |
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:
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.
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.
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.
Delphi has two native types of arrays:
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).
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;
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
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.
(in preparation)
(in preparation)
|
Scientific Programming |
©Nikolai V. Shokhirev, 2001-2005