|
Example - Interpolate
* The "\" symbol has been used at the end of lines to indicate that the line was too long for the display and has been
wrapped to the line below in order to fit within the width of the web page. These lines would need to be joined in the normal
program editor.
I. Introduction
This program will draw a smooth curve between points plotted on a graph by the user's mouse. Then, the area under the curve and length of the curve are calculated. The program demonstrates the use of matrices and vectors (for the interpolation) and integration (for calculating length and area).
II. Mathematics Involved
A curve is fitted to the points by using a series of quadratic curves of the form y=ax 2 +bx+c. We consider a set of three points at a time and fit a quadratic curve to these, then take the last of these three as the first of the next three points. However, to avoid a discontinuity in the overall curve's gradient where two sections meet, the derivative of the curves must be equal at this point as well as the position. We thus have the following equations for the points r to r+2:
where a' , b' and c' are the quadratic coefficients for the previous curve.
These can be solved uniquely for a , b and c by writing these formulae as a matrix equation:
 , where  .
Then we can invert the matrix to give  .
[The first values of a, b and c are determined by solving the equations
using a similar method, then the calculation proceeds iteratively.]
When all the equations for the lines have been calculated, the area under the overall curve is given by the sum of the definite integrals from beginning to end of each quadratic curve section.
The length of a line with equation y=f(x) and endpoints at x= xa and x= xb is given by  .
III. How the Program Works
1) Initialisation
A volatile file is created to hold the points. The file is indexed on the x field, the x-value of the point; this allows easy access to the points in order later on. The other fields hold the a, b and c values for the equation of the line, and the pixel position of the point. The beginning and endpoints of the line, (0,0) and (1,0) are stored in memory.
SUB createpointfile()
DIM x%
CREATE MEMORY "points"
ADD "x;num"
ADD "y;num"
ADD "a;num"
ADD "b;num"
ADD "c;num"
ADD "xpix;nmi"
ADD "ypix;nmi"
MAKE "points"
CREATE INDEX ON x.points
FOR x% = 0 TO 1 STEP 1
BLANK
x.points = x%
y.points = 0
a.points = 0
b.points = 0
c.points = 0
STORE
NEXT x%
END SUB
2) Adding of Points
Once the form has been set up, the main loop is called for adding and removing points. When the graph area is clicked with the left mouse button a point is added at that position. If the point is right-clicked then the point if it exists is removed from the graph.
IF ((foundit%% = 0) OR ( ABS (ypoints% - y%) > 0.004)) THEN
BLANK FILE "points"
x.points = x%
y.points = y%
xpix = mx%%
ypix = my%%
STORE
myform.DefFontBold = - 1
SET fc = myform.Add("","ellipse")
fc.Move(mx%% - 3,my%% - 3,6,6)
fc.OnRightClick = "DotRemove"
fc.tag = "dot"
END IF
FORM
END IF
3) The ProcessCurve Procedure
When the Complete Curve button is clicked, the ProcessCurve procedure is called. The following code shows the algorithm used to calculate the quadratic coefficients for the curves (the coefficients for the first curve having already been worked out).
WHILE (d%% < r%%)
' the first entry is the derivative at the first
' point of the curve
SELECT NEXT FILE "points" INDEX "x"
m.elements(1,1) = 2 * x.points
m.elements(1,2) = 1
m.elements(1,3) = 0
v.elements(1) = 2 * a% * x.points + b%
m.elements(2,1) = x.points * x.points
m.elements(2,2) = x.points
m.elements(2,3) = 1
v.elements(2) = y.points
SELECT NEXT FILE "points" INDEX "x"
m.elements(3,1) = x.points * x.points
m.elements(3,2) = x.points
m.elements(3,3) = 1
v.elements(3) = y.points
SELECT PREVIOUS FILE "points" INDEX "x"
v.multiply(m.invert())
a% = v.elements(1)
b% = v.elements(2)
c% = v.elements(3)
a.points = a%
b.points = b%
c.points = c%
STORE FILE "points"
d%% = d%% + 1
WEND
The curve is then drawn, and the integral calculations performed. For this we need an explicit equation for the line for every value of x between 0 and 1; the function func%() provides this. For any value of x, a record in the ‘points' file contains the quadratic coefficients that apply at that point. Using those coefficients an explicit value for the line function can be calculated for that x value.
FUNCTION func%(x%)
IF ((x% <= 0.0) OR (x% >= 1.0)) THEN
func% = 0.0
ELSE
' note that 'select key' should ensure
' that x.points >= x%
SELECT KEY x% FILE "points" INDEX "x"
IF (x.points > x%) THEN
SELECT PREVIOUS FILE "points" INDEX "x"
END IF
func% = (x% * x% * a.points) + (x% * b.points) + c.points
END IF
END FUNCTION
Similarly, the function funcdot% provides the value of the derivative of the line function at any point between 0 and 1.
FUNCTION funcdot%(x%)
IF ((x% <= 0.0) OR (x% >= 1.0)) THEN
funcdot% = 0.0
ELSE
' note that 'select key' should ensure
' that x.points >= x%
SELECT KEY x% FILE "points" INDEX "x"
IF (x.points > x%) THEN
SELECT PREVIOUS FILE "points" INDEX "x"
END IF
funcdot% = (2 * x% * a.points) + b.points
END IF
END FUNCTION
This enables the calculation of the area under the curve, and length:
' now do the area under the curve
' and the length of the curve
area$ = " ..."
length$ = " ..."
FORM
a% = integrate("func%(x%)","x%",0.0,1.0,0.0005,20)
area$ = STR$ (a%,"99.000")
FORM
b% = integrate("sqr(1 + (funcdot%(x%)^2))","x%",0.0001,0.9999,0.0005,20)
length$ = STR$ (b%,"99.000")
FORM
|
