[Back to GRAPHICS SWAG index] [Back to Main SWAG index] [Original]
Program Interpolation;
{
ÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛ
ÛÛÛÝÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÞÛÛÛ±±
ÛÛÛÝÛÛ ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ Least squares method interpolation ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ Aleksandar Dlabac ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ (C) 1995. Dlabac Bros. Company ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ ------------------------------ ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ adlabac@urcpg.urc.cg.ac.yu ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ adlabac@urcpg.pmf.cg.ac.yu ÛÛÞÛÛÛ±±
ÛÛÛÝÛÛ ÛÛÞÛÛÛ±±
ÛÛÛÝßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÞÛÛÛ±±
ÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛ±±
±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±
}
{
Interpolation is a technique used to calculate estimated function of
which we have some experimental points. Firstly, we have to choose order of
interpolation. Order of interpolation represents the order of estimated
function. For example, if order is 5 function will be in form of:
F(x) = a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0
This program (procedure Interpolate) returns coefficients a0..an, where
n is order. We have to determine the order considering positions of points.
}
Uses Crt, Graph;
Type Row = array [0..10] of real; { Maximal order }
Matrix = array [0..10] of Row; { is 10. }
Data = array [1..200] of integer; { Data for interpolation }
Coeff = array [0..10] of real; { Coefficients (also max. 10) }
Const Order = 5; { Order of interpolated curve. Should be smaller than }
{ number of points. For orders greater than 6-7 system }
{ (of equations) becomes unstable, so curve could be }
{ wrong. Also, for higher orders real type should be }
{ changed to extended, and $N switch should be on. }
Var Gd, Gm, I, J : integer;
R : real;
InversionOK : Boolean;
Mat : Matrix;
X, Y : Data;
C, YM : Coeff;
Function Power (Base,Pow:integer) : real;
Var Temp : real;
Begin
If (Base<0) and (Pow mod 2=1) then
Temp:=-Exp (Pow*Ln (-Base))
else
Temp:=Exp (Pow*Ln (Abs (Base)));
Power:=Temp
End;
Procedure MatrixInversion (Var A:Matrix; N:integer);
Var I, J, K : integer;
Factor : real;
Temp : Row;
B : Matrix;
Begin
InversionOK:=False;
For I:=0 to N do
For J:=0 to N do
If I=J then
B [I,J]:=1
else
B [I,J]:=0;
For I:=0 to N do
Begin
For J:=I+1 to N do
If Abs (A [I,I])<Abs (A [J,I]) then
Begin
Temp:=A [I];
A [I]:=A [J];
A [J]:=Temp;
Temp:=B [I];
B [I]:=B [J];
B [J]:=Temp
End;
If A [I,I]=0 then Exit;
Factor:=A [I,I];
For J:=N downto 0 do
Begin
B [I,J]:=B [I,J]/Factor;
A [I,J]:=A [I,J]/Factor
End;
For J:=I+1 to N do
Begin
Factor:=-A [J,I];
For K:=0 to N do
Begin
A [J,K]:=A [J,K]+A [I,K]*Factor;
B [J,K]:=B [J,K]+B [I,K]*Factor
End
End
End;
For I:=N downto 1 do
Begin
For J:=I-1 downto 0 do
Begin
Factor:=-A [J,I];
For K:=0 to N do
Begin
A [J,K]:=A [J,K]+A [I,K]*Factor;
B [J,K]:=B [J,K]+B [I,K]*Factor
End
End
End;
A:=B;
InversionOK:=True
End;
Procedure Interpolate (X,Y:Data; var A:Coeff);
Var I, J, K : integer;
Begin
For I:=0 to Order do
For J:=0 to Order do
Mat [I,J]:=0;
For I:=0 to Order do
For J:=0 to Order do
For K:=1 to 10 do
Mat [I,J]:=Mat [I,J]+Power (X [K],I+J);
MatrixInversion (Mat,Order);
For I:=0 to Order do
Begin
YM [I]:=0;
A [I]:=0
End;
For I:=0 to Order do
For J:=1 to 10 do
YM [I]:=YM [I]+Y [J]*Power (X [J],I);
For I:=0 to Order do
For J:=0 to Order do
A [I]:=A [I]+Mat [I,J]*YM [J];
End;
Begin
DetectGraph (Gd,Gm);
InitGraph (Gd,Gm,'');
X [1]:=100; Y [1]:=100;
X [2]:=150; Y [2]:=180;
X [3]:=200; Y [3]:=250;
X [4]:=250; Y [4]:=270;
X [5]:=300; Y [5]:=300;
X [6]:=350; Y [6]:=320;
X [7]:=400; Y [7]:=290;
X [8]:=450; Y [8]:=230;
X [9]:=500; Y [9]:=130;
X [10]:=550; Y [10]:=90;
Interpolate (X,Y,C);
SetFillStyle (SolidFill,Red);
For I:=1 to 10 do {< Draws red points (input }
Bar (X [I]-1,Y [I]-1,X [I]+1,Y [I]+1); {< values - X & Y). }
SetColor (Green);
For I:=10 to 55 do {< Draws calculated curve. }
Begin
R:=0;
For J:=0 to Order do
R:=R+C [J]*Power (I*10,J);
If I=10 then
MoveTo (I*10,Round (R))
else
LineTo (I*10,Round (R))
End;
Repeat Until ReadKey<>#0;
CloseGraph
End.
[Back to GRAPHICS SWAG index] [Back to Main SWAG index] [Original]