The magazine of the Melbourne PC User Group

SOLV: A Solver of Equations in One Unknown
Max G. Chandivert

In this article I would like to use a small program to illustrate four propositions:
  • The power in your machine is available to you without having to use a large and expensive commercial program.
  • Writing a small program for a particular purpose is simple, inexpensive and challenging.
  • Such a small program can be written in interpreted BASIC (GWBASIC or BASICA) at NO cost, except for some programming time.
  • The BASIC interpreter is, in some respects, more versatile than a compiler as it can accommodate self-modifying code. As an example let us look at the problem of solving an equation in one unknown. Surely the machine for which yourself or your firm has just paid some $2000 can do what some pocket calculators do routinely through their "SOLVE" key!
Problem 

We will assume that the expression in which you are interested has been put into the form Y = F(X), and that you are seeking the roots of Y, i.e. the value(s) of X which make(s) Y vanish. The SECANT formula is the basis for the method presented here. We start with two values of X, A and B, supplied by the user, and calculate the corresponding values of Y, YA and YB. The user may use any information he/she has about Y to pick A and B so that YA and YB are of opposite signs, but this is not critical. Then in the absence of discontinuities at least one root must be located between A & B. We calculate the slope of the secant joining points (A, YA) and (B, YB), and mathematically seek its inter-section with the X axis. That new guess X = G yields anew value of Y, YG, which is hopefully closer to the root required. Then we substitute B to A, G to B and go around the loop again. A limit is placed on the number of iterations so as not to waste time if the process does not converge quickly, which is rare. 

Naturally not always do the initial guesses A, B yield values of Y which are +, -. We might have -, - or +, +. This is where some "art" intrudes into an otherwise scientific subject. In particular when any two values YA, YB have close magnitudes, the joining secant would be nearly parallel to the X axis and prevent convergence of the process to a root. This situation must be detected by the program and remedied by bending the secant appropriately. 

A particularly difficult case is that of an equation having a double root at one point, i.e. where the graph of Y is tangent to the X axis, and in the area lies wholly on one side of it. In this case the program introduces the signs of the derivatives of Y at (A, YA) and (B, YB) to control the program logic and substitutes an arbitrary guess to that calculated by the secant formula. Naturally the code to cater for those "special" conditions has to be merged to that catering to the "normal" condition with as little interference as possible.

Solution 

The listing SOLVBAS is a small program which is surprisingly versatile in finding roots of equations in one unknown. Some examples of equations are given below with suggestions for trial values. To run the program, load it under the BASIC interpreter, GWBASIC or BASICA, and press F2; the rest is self explanatory. The program creates a small temporary file in the current directory, so it is necessary to have a free cluster there. After each solving attempt the user is asked whether he/she wishes to repeat with different initial guesses or to change the equation to be solved. The maximum number of iterations is hard-coded to 50 in the present version.

 SOLV.BAS  [Source listing may be downloaded from this link] 

100 '- SOLV.BAS-A program to solve equations in one unknown -
110 KEY OFF : SCREEN 0: COLOR 14, 1 : CLS: DEFINT I-N
120 LOCATE 2, 20: PRINT "- SOLV - SOLVER OF EQUATIONS IN ONE UNKNOWN -":PRINT
130 PRINT SPC(12);"by Max Chandivert - MS/PC-DOS Interpreter version - 900822"
140 A$ ="6*X^5 -13*X^4 +27*X*X*X +7*X*X -29*X +10"
150 PRINT: PRINT "F(X) = "A$
160 PRINT "Is another expression required Y/N? : .: R$ = INPUT$(1)
170 PRINT R$: IF R$ = "N" OR R$ = "n" THEN 220
180 INPUT "New expression in BASIC code: F(X) _ ",A$
190 F$ = "TTT.BAS": OPEN F$ FOR OUTPUT AS #1
200 PRINT #1, "140 A$ ="+CHR$(34)+A$+CHR$(34): CLOSE #1
210 CHAIN MERGE F$, 140
220 F$ = "TTT.BAS"
230 OPEN F$ FOR OUTPUT AS #1: PRINT #l, "250 DEF FNF(X)="+A$: CLOSE #1
240 CHAIN MERGE F$, 250
250 DEF FNF(X)=6*X^5 -13*X^4 +27*X*X*X +7*X*X -29*X +10
260 TOL1=.000001: TOL2=1!: PRINT: PRINT "TOL1, TOL2: ";TOL1,TOL2
270 PRINT "Do you wish to change those values, Y/N? : ,: R$ = INPUT$(1)
280 PRINT R$
290 IF R$ _ "Y" OR R$ = "y" THEN INPUT "Enter new values: ",TOL1, TOL2
300 PRINT: INPUT "Enter A,B trial values for X: ",A,B
310 IF A<B THEN TEMP=A: A=B: B=TEMP
320 YA = FNF(A): YB = FNF(B): ITAG = 0
330 PRINT :PRINT "Initial values of Y are YA,YB= "YA, YB
340 IF ABS(YA) < TOLl THEN G = A: N = 0: GOTO 540
350 IF ABS(YB) < TOLI THEN G = B: N = 0: GOTO 540
360 PRINT : PRINT "Successive trial values of variable, and resulting Y."
370 FOR N=1 TO 50
380 I=SGN(YA)*SGN(YB): SLOPE = (YB-YA)/(B-A): G = B-YB/SLOPE
390 YG = FNF(G) ': PRINT N, I, ITAG To facilitate development
400 IF ABS(YG) < TOLl THEN 540
410 IF ABS(YG) < ABS(YB) THEN 510
420 IF I=1 THEN ITAG=ITAG+l ELSE ITAG=0: GOTO 510 'YA, YB are opposite  signs
430 IF ITAG > 2 THEN 450 ' Now YA, YB have same sign
440 IF ABS(SLOPE)<TOL2 THEN G = B - YB/TOL2: GOTO 490 ELSE GOTO 510
450 YAA = FNF(A+.001*ABS(A)): YBB = FNF(B+.001*ABS(B))
460 IF YAA>YA THEN IDA =1 ELSE IDA --1
470 IF YBB>YB THEN IDB -1 ELSE IDB --1
480 IF IDA*IDB = -1 THEN G = .5*(A+B) ELSE 510
490 YG = FNF(G)
500 IF ABS(YG)<TOLl THEN 540
510 A=B: B=G: YA=YB: YB=YG
520 PRINT G, YG: NEXT
530 PRINT "N= "N"-NO SOLUTION FOUND": GOTO 560
540 PRINT "Y within ";TOL1:" of ZERO for X= ";G
550 PRINT "There was ";N;" iteration(s)."
560 PRINT: PRINT "Try other starting values A, B, Y/N? : ";: R$ =  INPUT$(1)
570 PRINT R$: IF R$ = "Y" OR R$ = "y" THEN 300
580 PRINT "Solve another equation, Y/N? : " ;: R$ = INPUT'$(1): PRINT R$
590 IF R$ = "Y" OR R$ = "y" THEN CLS: GOTO "183 ELSE KILL "TTT.BAS": END

  Reprinted from the October 1990 issue of PC Update, the magazine of Melbourne PC User Group, Australia