home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turbo Toolbox
/
Turbo_Toolbox.iso
/
sonderh1
/
fitpoly.pas
< prev
next >
Wrap
Pascal/Delphi Source File
|
1986-11-25
|
7KB
|
266 lines
(*-------------------------------------------------------------------------*)
(* Programm zur Berechnung von Ausgleichspolynomen.
Die x/y-Wertepaare, zu denen die Polynome berechnet werden sollen,
werden aus einem Textfile eingelesen.
Zum Verfahren s.:
Stoer J. , Heidelberger Taschenbuecher Nr.105:
" Einfuehrung in die numerische Mathematik I ",
S. 168 ff,
Springer-Verlag , Berlin...
Das Programm wurde auf einem ATARI 520ST+ mit
Pascal ST plus von CCD geschrieben. *)
(*-------------------------------------------------------------------------*)
PROGRAM poly_fit(input,output);
CONST
maxzeil = 200;
maxspalt= 20;
TYPE
(* Manche Compiler koennen ueber dyn. Var. groessere Felder erzeugen.
Der von CCD allerdings nicht. *)
ahi = ARRAY [1..maxzeil,1..maxspalt] of REAL;
ptr = ^ahi;
VAR
zeil, (* Speichert die Anzahl der Gleichungen *)
spalt, (* " " " " Unbekannten *)
grad : INTEGER; (* " den Grad des Fitpolynoms *)
a : ptr; (* Feld fuer die Koeffizienten des Gleichungssystems *)
d, (* Hilfsfeld *)
l : ARRAY [1..maxspalt] of REAL; (* Feld fuer das Ergebnis der Rechnung *)
singmat : BOOLEAN; (* Matrix singulaer oder nicht ?! *)
(*-------------------------------------------------------------------------*)
PROCEDURE householder;
VAR
sigma, s, sum, beta : REAL;
i, j, k : INTEGER;
BEGIN
Write ('Rechne');
singmat := FALSE;
j := 1;
WHILE (j <= spalt) AND NOT(singmat) DO
BEGIN
Write ('.'); (* Damit man sieht, dass etwas passiert *)
sigma := 0;
FOR i := j TO zeil DO
sigma := sigma + a^[i,j] * a^[i,j];
IF sigma <> 0 THEN
BEGIN
IF a^[j,j] < 0 THEN
BEGIN
s := Sqrt (sigma); d[j]:=s;
END
ELSE
BEGIN
s := -Sqrt (sigma); d[j]:=s;
END;
beta := 1/(s*a^[j,j]-sigma); a^[j,j] := a^[j,j] - s;
FOR k := j+1 TO spalt+1 DO
BEGIN
sum := 0;
FOR i := j TO zeil DO
sum := sum + a^[i,j] * a^[i,k];
sum := beta * sum ;
FOR i := j TO zeil DO
a^[i,k] := a^[i,k] + a^[i,j] * sum;
END;
END
ELSE
BEGIN
WriteLn; WriteLn;
WriteLn ('Singulaere Matrix');
singmat := TRUE;
END;
j := Succ (j);
END;
WriteLn;
END; (* householder *)
(* ------------------------------------------------------------------------- *)
PROCEDURE ruecktransform;
VAR
i, j : INTEGER;
BEGIN
FOR i := spalt DOWNTO 1 DO
BEGIN
l[i] := a^[i,spalt+1];
FOR j := spalt DOWNTO i+1 DO
BEGIN
l[i] := l[i] - a^[i,j] * l[j];
END;
l[i] := l[i]/d[i];
END;
END; (* ruecktransform *)
(* ------------------------------------------------------------------------- *)
PROCEDURE pagen; (* Bildschirm loeschen und Cursor nach links oben *)
BEGIN
Write (Chr(27));
Write ('E');
END; (* pagen *)
(* ------------------------------------------------------------------------- *)
PROCEDURE ende;
VAR
anything : CHAR;
BEGIN
Write ('Mit <RETURN> zurueck zum Desktop');
Read (anything);
END; (* ende *)
(* ------------------------------------------------------------------------- *)
PROCEDURE ausgabe;
VAR
i : INTEGER;
BEGIN
WriteLn;
IF zeil > spalt THEN
Write('Koeffizienten des Fitpolynoms: ');
FOR i := 1 TO spalt DO
BEGIN
IF (i-1) MOD 3 = 0 THEN
BEGIN
WriteLn;
END;
Write('a[',Pred(i):2,']',' = ', l[i]:14,' ');
END;
WriteLn;
WriteLn;
END;
(*-------------------------------------------------------------------------*)
(* Diese Prozedur kann auch so geschrieben werden, dass kein
zweimaliges Lesen des Datenfiles noetig ist. Das wird dann aber
eher unuebersichtlicher. *)
PROCEDURE einlesen;
VAR
name : STRING [25];
tex : text;
n,i,j,numcount : INTEGER;
buff : STRING [240];
c1,c2 : CHAR;
BEGIN
pagen;
Write('Name des Datenfiles: ');
ReadLn(name);
Write('Grad des Fitpolynoms: ');
ReadLn(grad);
WriteLn;
Write('Lese Daten...');
Reset(tex,name);
REPEAT (* BEGIN suchen *)
ReadLn(tex,buff)
UNTIL pos('BEGIN',buff)=1;
numcount:=0;
REPEAT
ReadLn(tex,buff);
numcount := Succ(numcount); (* Zeilen zaehlen *)
UNTIL Pos('END',buff) = 1;
Reset(tex,name); (* Zurueck an den Anfang *)
REPEAT
ReadLn(tex,buff)
UNTIL Pos('BEGIN',buff) = 1;
i := 0;
n := 0;
FOR n := 1 TO Pred(numcount) DO
BEGIN (* Koeffizienten einlesen *)
a^[n,1] := 1; (* erste Spalte ueberall 1 *)
Read(tex,a^[n,2]); (* x-Wert einlesen und Potenzen *)
i := 3;
WHILE i<= Succ(grad) DO
BEGIN (* ... berechnen. *)
a^[n,i] := a^[n,2] * a^[n,Pred(i)];
i := Succ(i);
END;
Read(tex,a^[n,grad+2]);
Write('.');
END;
zeil := numcount; (* Zeilenzahl festlegen *)
spalt := Succ(grad); (* Spaltenzahl " *)
WriteLn;
WriteLn;
END;
(*-------------------------------------------------------------------------*)
(* Diese Prozedur schreibt das Fitpolynom als komplette Pascal-
Function in einen Textfile. Von dort kann es dann in eigenen
Programme eingebaut werden. *)
PROCEDURE polynom_wegschreiben;
VAR
ch : CHAR;
name : STRING [25];
tex : TEXT;
i : INTEGER;
BEGIN
REPEAT
WriteLn
('Soll das Polynom in eine Textdatei als Pascalfunktion geschrieben werden?');
Write('(J/N)');
Read(ch);
WriteLn;
UNTIL ch IN ['j','J','n','N'];
IF ch IN ['j','J'] THEN
BEGIN
Write('Name der Ausgabedatei? ');
ReadLn(name);
ReWrite(tex,name);
WriteLn(tex,'function ywert(x : REAL) : REAL;');
WriteLn(tex);
WriteLn(tex,'VAR');
WriteLn(tex,' y : REAL;');
WriteLn(tex);
WriteLn(tex,'BEGIN');
WriteLn(tex,' y := ',l[Succ(grad)]:15:12,';' );
FOR i:=grad DOWNTO 1 DO
BEGIN
Write(tex,' y := y * x ');
IF l[i]>0 THEN
Write(tex,'+ ')
ELSE
Write(tex,'- ');
WriteLn(tex,ABS(l[i]):15:12,';');
END;
WriteLn(tex,' ywert := y;');
WriteLn(tex,'END;');
END;
END;
(*-------------------------------------------------------------------------*)
BEGIN (* poly_fit *)
New(a);
einlesen;
householder;
IF NOT(singmat) THEN
BEGIN
ruecktransform;
ausgabe;
polynom_wegschreiben;
END;
ende;
END.