Math Development Program, Part 4

With the continuing upheaval in the computer industry, one can only speculate half-heartedly about the future of computers. They won’t go away, but how many of the different types will survive? How advanced will they get and how much will they cost? Whatever those answers may be, the reasons for using computers remain the same. Two of the current philosophies are string manipulation (word processing, communications) and number manipulation (spreadsheets, games). The computer thrashes blindly through memory locations doing its’ job. The difference between two programs lies in the sequence of execution or algorithm, if you will. The crux here is, how do you write the perfect algorithm for your given job? First, you need a job to do. You then must realize that a computer is just a tool; it only does what you tell it and finally; practice, practice, practice!!

Knowing how to program does not get you very far any more. My five year old son can do very basic programming. You need more, much more. A solid math background is a good start regardless of other expertise. A good graphic game that has more than straight line movement requires some premeditated algorithm to accomplish that move. Once the computation is made, the move can be stored and regenerated instantly with minimum effort. If you don’t know the algorithm then forget it. Your game’s a bore. Computer simulations are pre-planned, they don’t just happen. People slave for months (right Fred?) over a great program.

Dissemination of information is what this magazine is all about. Many computer magazines support nothing but product reviews-all bun, no beef! In my series on math, some of the routines were translated directly from Fortran (Gauss Elim.) while others were original (Cramers rule any size, interpolation) along with plenty of space saving tips. SE in its’ glory is a 20K ASCII file in 12.5K Sinclair basic (12 page listing). There is still plenty to be done. Still, no matter how many subroutines you have at your disposal, you need some knowledge of how to use them. My philosophy all along has been to explain the basics of the algorithms in the programs. They were all hand derived at some point in the past. The last major math installment follows; Polynomial and real time differentiation (derivative).

Differentiation

This is a vastly important subject. Your mind is extremely fast at calculating a derivative. Every time you catch a ball or stop for a traffic light, your mind is busy differentiating, among other things. Differentiation is the opposite of Integration. An integral covers an area or is total work as noted by addition (miles, minutes, your electric meter reading) or multiplication (eg. sq. ft., ft. lbs., mile-minutes, volt-amps, kilowatthours) where differentiation is denoted by division (eg., velocity or speed in ft/sec, or force in lbs./inch). To make a distinction, acceleration is a second derivative (ft/sec/sec), or the derivative of the first derivative (velocity). In other words, speed changes with time. Totally confused yet? Sorry, just note that derivatives are a small piece of what is happening to the whole function at a given moment.

When you have a lot of data, it is easy to plot and see what is going on with the trends and tendencies of your information. If you have precious few bits of information, then it is very helpful to know a lot about those few points. The derivative is the means by which you develop information about your points. In simple terms, on an X,Y coordinate system, the derivative is the rise over the run (Y/X). It is noted as dy/dx or delta Y/delta X. These two terms have different meanings, but are identical for linear models. These delta and dx functions stand specifically for the change in Y with respect to X. X will usually be time for real time data, but it can be anything we want it to be. Time is nice since it is impossible to sample two data points and have no dx (change in time).

Should you get out your old (or new) Calculus book, you can look up the formula for calculating a derivative. It looks something like this:

f'(X)=((f(X+dx)-f(X))/dx)

You may find this a little easier to see:

dy=(Y(I)-Y(I-1))/(X(I)-X(I-1))

If you are familiar with trigonometry, then you may recall that Y/X or dy/dx will give the tangent of the angle between the X axis and a straight line denoting the change in Y values. This makes for easy plotting of data by summing values (Logo like). Since Y values are plotted directly above the X values then,

Y=Y(I-1)+(X(I)-X(I-1))*TAN((Y(I-1)/(X(I)-X(I-1))

Now, you may think that this is just a bit cumbersome to calculate for each value you read (that’s INPUT in Sinclairese) in, but after the first time through it becomes much easier. Y becomes Y(I-1) and X(I)-X(I-1) can be dx. The above formula simplifies to:

Y=Y+dx*TAN(Y/dx)

You then need only PLOT X,Y (provided they are on scale). This method is still somewhat inconvenient and all it really does (in this form) is guess at what the new Y value is. The next Y value to be read would be the actual Y value. If these values are not equal, then you are not dealing with a linear function. This does shed some light on whether your function is increasing at a faster rate or a slower rate. This is in fact some indication of the second derivative. An example would be a time, distance vector. If time values are plotted on the X axis and Y values are the distance you travel after a certain time, then it is easy to see that at a constant speed (a linear-straight line function) you can predict very easily where you will be at any time. Speed is the first derivative. At 60 (MPH), you will move 88 feet in the next second. Air Force bombardiers thrive on first derivatives. However, if you have a lead foot like myself and accelerate until it is necessary to brake to avoid hitting the car in front, then you can not accurately predict distances versus time. Bombardiers aren’t thrilled about second derivatives.

One more thing to confuse you about these functions is the sign of each (first and second) type of derivative is independent of the other. If both derivatives are positive, then the Y values are getting larger faster. If the first is positive and the second is negative, then the Y values are getting larger slower. If the first is positive and the second is zero, then all Y values in that region should lie on a straight line (increasing distance-slope up). Similarly, if both are negative, then you are losing ground faster and faster. With the second derivative positive, you are slowing down (but still going backwards). Again, if the second derivative is zero while the first is negative, the speed is constant (linear, backwards or decreasing-slope down). If the first derivative is zero, this implies the the distance is constant or, you have stopped.

Nothing to it, right? Any good math book will say that that is just how derivatives operate. The engineer will say that it sure looks good on paper. Don’t program your robot to pour a cup of coffee this way though. The problem remains that the source of the data is subject to spikes and transmission errors. Even though they may be one in a billion, bad data can make that coffee fall right in your lap! Taking two points to evaluate a derivative implies a mid-point value and not a true real time value. The solution to high speed related rates type problems lies in not only software, but hardware as well. Fast data aquisition may be vital. However, if your robot is not running a horse race, then software can do the job. Taking two data points very close together may get you closer to a real time derivative but it will not improve your accuracy too much and you are still subject too data error. Taking more data points will help much more although your true evaluation point may be slightly further from a present time position.

Why worry at all about the derivative? What practical use is it? The problem of starting, moving and stopping a robot smoothly has plagued AI (artificial intelligence) scientists for years. There are many sophisticated schemes used to manipulate robots. Still, the manual ‘LEARN’ mode prevails at this time. Predetermined algorithms are being developed right now. The best AI books available are just collections of papers by various experts. (This is still a wide open field for anyone to partake). Getting back to basics, data acquisition and display can be handled conveniently with present day knowledge. Two easy ways to keep track of relative changes in your data are moving trend line (eg., stock market) and polynomial differentiation. There are two methods of differentiating polynomials as well. One is with a few points and the other is with a function (of course). If data input is sporadic, then a trend line is advantageous. On the other hand, regular data input (without lots of spikes and noise) can take advantage of an N-Point (differentiating an interpolating polynomial) algorithm. If you desire to study the data by collecting and ‘Post-Processing’, then fitting and differentiating the function will do nicely.

The first table is composed of the N-Point functions. The more points, the more accuracy. Three and five points are particularly nice because the mid-point (which happens to be the point you are evaluating the derivative for) drops out and leaves one less calculation. Besides, 4 point actually require 7 values although three drop out of the equation. The convention used in the following table has (I) as the point for which the derivative is being calculated for. Remember that the values are interpolated. Evaluating trigonometric functions like sines and cosines requires more than 4 points. In SWN 1/3, the error for interpolating holds for differentiating also. For N=4 and N=5 points, there are third derivatives and N=5 has a fourth derivative. If you have interest in these, send us a SASE and we’ll send them to you.

N-POINT DIFFERENTIATION
IN THE MIDDLE OF AN INTERVAL

N=3   dy=(-Y(I-1)+Y(I+1))/(2*dx)
      d2y=(Y(I-1)-(2*Y(I))+Y(I+1))/(dx**2)

N=4   dy=(Y(I-3)-27*Y(I-1)+27*Y(I+1)-Y(I+3))/(48*dx)
      d2y=(Y(I-3)-Y(I-1)-Y(I+1)+Y(I+3))/(2*dx**2)

N=5   dy=(Y(I-2)-8*Y(I-1)+8*Y(I+1)-Y(I+2))/(12*dx)
      d2y=(-Y(I-2)+16*Y(I-1)-30*Y(I)+16*Y(I+1)-Y(I+2))/(12*dx**2)

As in past issues, I always like to give more than one way to accomplish a task if possible. Different circumstances require a change of tactics. Interpolating, means estimating an exact course between known values. These estimations can (but not necessarily) gyrate strangely. As in least squares fitting (SWN 1/3), if you are interested in the trend of the change over a length of time or distance, a simple least squares fit will do. When I say simple, I mean just that. If data is coming in at set intervals, then dx (as well as all X) is a constant. Once the interval is defined, over half of the least square calculations drop out. It becomes necessary to keep track of only Y values. If you only want to keep track of a given amount of data and can discard the rest, then a moving trend line will give a very stable derivative (the stock market should be so lucky). To accomplish this, just divide or subtract out the old Y value and multiply or add in the new one. Recalculate the slope (M), the intercept (B) will become a pivot (constant), and you have an indication of the first derivative. The values will be relative but will show the trend nicely. If you save the old M, then subtracting the old M from the new M will give you an indication of the second derivative. Simple and fast! One hitch, after subtracting, dividing, etc., on the same set of numbers for long periods of time, watch out for roundoff error. Every so often you should reinitialize your array and start fresh once again. If your X values are at varying times, then you must redo the whole least squares after each point. Cramer’s Rule is pretty quick though. The listing given runs in about 0.5 seconds in fast, 3.6 in slow. That includes those time consuming extra scrolls and prints also. This technique is certainly not instantaneous, but spikes in the incoming data will be averaged out quite nicely. If the data rate is fairly slow, then higher orders of derivation (a bigger set of simultaneous equations – remember them?) can be employed to get a more accurate derivative. Least squares routines utilize the partial derivatives of each degree of a function to evaluate a fit. Any degree of fit can be acheived if wanted or needed. This method only gives the trend through the entire fit area. However, if the function is developed, then the next little technique will give the instantaneous derivative for any spot inbetween the first and last points defined.

Horner’s Rule

To ease the implementation of this technique, we must set up our function so as to evaluate it with the utmost speed and efficiency. The technique is called Horner’s Rule (HR). This simple method reduces polynomial calculations by a third and gets rid of those slow exponent evaluations. We can set up a loop or fill a string with a function. The loop method is more dynamic, but the string method will retain the function for whatever else you can cook up.

Use Horner’s Rule wherever it fits; it is just good programming.

This very short routine takes the coefficient (X(I))of the highest power times the value of the point you are evaluating, adds it to the next lower power X(I) and loops back to do it again. This continues for all terms.

The Horners Rule conversion looks like this:

4*X**3+2*X**2-3*X+1 becomes the following,

((4*X+2)*X-3)*X+1

The first time through the loop in line 9550, will give the Y value or is the same as using VAL F$. I believe that this routine is faster than VAL since there are no powers, but I’m not sure. Besides, the following times through the loop gives all of the derivatives. The swap terms routine is necessary to get the powers lined up in the right order to take advantage of the routine. If you use the interpolation routine, then you can bypass the swap routine completely. The create function routine is just to compare the Val function with the HR routine. The highest power coefficient should be entered as X(1) and so on, and they must be present before you ‘create’ your function. A function in F$ is not necessary to evaluate the derivative though. Just put X in an incrementinng loop and GOSUB listing 3. you can watch all the derivatives change with each point evaluated.

I certainly hope you have enjoyed and comprehended this Math series. I know how difficult this material is to learn. Believe me, it takes a lot of patience and PRACTICE!

Differentiation Demo

10 REM SINE TREND
20 SAVE "MATH[illegible]"
30 REM LISTING NO. 2
40 LET N=4
50 DIM Y(N)
60 LET M=0
70 SLOW
90 REM TIME INTERVAL
100 LET T=3
110 LET T1=0
120 LET T2=0
130 LET Y=0
140 LET TY=0
150 LET X=0
160 FOR I=1 TO N
170 LET T1=T1+T
180 LET T2=T2+T*T
190 LET Y(I)=0
200 LET Y=Y+Y(I)
210 LET TY=TY+T*Y(I)
220 LET T=T+T
230 NEXT I
240 LET CR=N*T2-T1*T1
250 IF NOT CR THEN STOP
260 REM ROUTINE STARTS HERE
270 FOR J=1 TO N
280 LET M1=M
290 LET M=(N*TY-Y*T1)/CR
300 LET M2=M-M1
310 LET Y=Y-Y(J)
320 LET TY=TY-T*Y(J)
330 LET X=X+3
340 SCROLL
350 SCROLL
360 REM PRINT OR PLOT Y,M AND M2
370 PRINT Y(J),M
380 SCROLL
390 SCROLL
400 PRINT ,M2
410 REM LET Y(J)=USR INPUT (Y(J)) VAL
420 LET Y(J)=SIN (X/180*PI)
430 LET Y=Y+Y(J)
440 LET TY=TY+T*Y(J)
450 REM PEEK FRAMES IN SLOW TO GET LOOP TIME (SECONDS)
460 REM PRINT (PEEK 16436+256*PEEK 16437)/-60+1092.25
470 REM MULTIPLY BY 6.75 FOR FAST MODE
480 REM PAUSE DELAY TIME OR GOSUB OTHER ROUTINE
490 NEXT J
500 GOTO 270
510 REM
520 REM LISTING NO. 1
530 REM REAL TIME DIFFERENTIATION
540 REM LET DX=1
550 LET DX=30*PI/180
560 SLOW
570 FOR I=0 TO 90 STEP 10
580 LET X=I*PI/180
590 LET DY1=(SIN (X-2*DX)-8*SIN (X-DX)+8*SIN (X+DX)-SIN (X+2*DX))/(12*DX)
600 LET DY2=(-SIN (X-2*DX)+16*SIN (X-DX)-30*SIN X+16*SIN (X+DX)-SIN (X+2*DX))/(12*DX*DX)
610 LET ER1=DY1-COS X
620 LET ER2=DY2-(-SIN X)
630 SCROLL
640 PRINT "DEG  1ST. DER."," ERROR"
650 SCROLL
660 PRINT I;"    ";DY1
670 SCROLL
680 PRINT ,ER1
690 SCROLL
700 PRINT "2ND.  DER."," ERROR"
710 SCROLL
720 PRINT DY2,ER2
730 SCROLL
740 NEXT I
750 STOP
760 REM
4740 REM LISTING NO. 3
4748 REM FROM SE
4749 REM (1 DIM )SWAP TERMS
4750 LET L=M
4760 FOR I=A TO M/2
4770 LET T=X(I)
4780 LET X(I)=X(L)
4790 LET X(L)=T
4800 LET L=L-1
4810 NEXT I
4814 REM INSERT NEXT LINE IN SE
4815 IF Z=7 THEN RETURN
4819 REM CREATE A POLYNOMIAL FUNCTION (FROM SE)
4820 LET F$=""
4830 FOR I=A TO M
4850 LET F$=F$+STR$ X(I)+"*X**"+STR$ (I-1)+" "
4860 IF I<>M THEN LET F$=F$+"+"
4870 NEXT I
4879 REM CHANGE 4880 TO :
4880 IF NOT AL THEN GOTO (5486 AND Z<7)+(9260 AND Z=7)
8998 REM
8999 REM
9000 REM POLYNOMIAL DERIVATIVE
9100 REM INPUT F$ OR FIT WITH SE PROGRAM
9110 REM INPUT COEFFICIENTS AS X(J) OR GET FROM SE
9230 LET Z=7
9240 LET AL=NOT PI
9250 REM GOTO 4820
9260 REM HORNERS RULE IN A LOOP
9265 PRINT F$
9270 DIM B(M)
9280 GOSUB 4750
9350 REM DERIVATIVES
9500 FOR I=1 TO M
9510 LET B(I)=X(I)
9520 NEXT I
9530 FOR J=1 TO M
9540 FOR I=2 TO M+1-J
9550 LET B(I)=B(I)+B(I-1)*X
9560 NEXT I
9570 PRINT "DER. ";J-1;" AT X = ";X;" = ";B(I-1)
9580 NEXT J
9590 GOSUB 4750
9600 RETURN