This is the source code which is discussed in the ITERATION PAPER.
It is written for Microsoft's QB45 syntax and can be used as is within that language environment.
Follow these steps,
1. Save this webpage using your browser's functions.
2. Use a plain ASCII text editor, WINDOW's Notepad or DOS Editor, and delete the HTML code lines along with these instructions.
There are three lines at the end of the file which should also be deleted.
3. Save the modified file as "ITERCODE.BAS", or other name of your choice, in your QB45 directory or whatever other basic language environment you use.
' Delete the above lines and this one.
' RE-SAVE "ITERCODE.BAS"
' Last STORED on Jul 21, 97 at 3:30
' Copyright by Jerry Hughes 1997
' Modification history:
' 1997
' Mar 24: This program accompanies the ITERATION paper which describes the
' iteration principles used in this code.
' The main program, including some GOSUB routines, contains code
' to aid in data entry and change as well as solution results.
' The structure of the code employs FUNCTION definitions where the
' user is allowed to define application equations and code.
' The syntax of this code is for MicroSoft's QBASIC 4.5 and use
' of their QBASIC software is required to use this code as is. Otherwise,
' translate the code into your basic's syntax or programming language.
'
' A full and complete DISCLAIMER is declared by the composer.
' Use at your own risk; verify results and applicablity.
'
DEFDBL X-Z
DECLARE FUNCTION FUNC1# (x, D$)
DECLARE FUNCTION FUNC2# (x, D$)
DECLARE FUNCTION FUNC3# (x, D$)
DECLARE FUNCTION FUNC4# (x, D$)
DECLARE FUNCTION FUNC5# (x, D$)
DECLARE FUNCTION LGT (x)
DIM xOFy(10), Y1eqY2atX(10), y1atx(10)
DIM LYmn(10), LXmn(10), LYmx(10), LXmx(10), FPy(10), FPx(10)
It0$ = "0": It1$ = "1": It2$ = "2": It3$ = "3": It4$ = "4": It5$ = "5": It6$ = "6"
It7$ = "7": It8$ = "8": It9$ = "9": ItA$ = "A": ItB$ = "B": ItC$ = "C"
ItD$ = "D": ItE$ = "E": ItF$ = "F": ItG$ = "G": ItH$ = "H": ItI$ = "I"
ItJ$ = "J": ItK$ = "K": ItL$ = "L": ItM$ = "M": ItN$ = "N": ItO$ = "O"
Ast$ = "*": Null$ = "": No$ = "No ": Yes$ = "Yes"
NotDone$ = "Not Done": Done$ = "Done"
Anal$(1) = "Inverse ": Anal$(2) = "Mx/Mn/FP ": Anal$(3) = "Intersect "
AnaDone$ = NotDone$
Fmt00$ = "##.####^^^^"
Fmt01$ = "& ##.####^^^^"
FOR I = 1 TO 79: Line$ = Line$ + " ": NEXT I
BlkHL$ = MID$(Line$, 1, 39)
Blu = 1: Gre = 2: Cya = 3: Red = 4: Mag = 5: Bro = 6: Whi = 7: Gra = 8
Bblu = 9: Bgre = 10: Bcya = 11: Bred = 12: Bmag = 13: Byel = 14: Bgra = 15
LColor = 9
Mx = 21: Ip = 22: BotLin = 23 ' Note/message, Prompt, Bottom lines
PI = ATN(1) * 4
Rpd = PI / 180
' **********************************
' Default values for input variables
' Define these values for intitial startup
' **********************************
FeqMax = 5 ' Maximum # of FUNCTIONs defined
Feq0 = 2: x = .000125: Feq = Feq0 ' Function to analize
GOSUB Selfunc: FDesc0$ = FDesc$ ' Function description
DoMin = 0: DoMax = 7 ' Domain
Npts = 100 ' Number of increments in domain
Anal = 1 ' Type of analysis, inverse
YofX = 10# ' Given function value
Tol# = .000001 ' Tolerance
AborTime = 90 ' Time out limit during iterations
Feq1 = 2: x = .000125: Feq = Feq1 ' First function for intersection
GOSUB Selfunc: FDesc1$ = FDesc$ ' and its description
Feq2 = 3: x = .000125: Feq = Feq2 ' Second function for intersection
GOSUB Selfunc: FDesc2$ = FDesc$ ' and its description
Progress = 1 ' Show progress status
ProgStat$ = "On"
' ************************************
' ************************************
GOTO Options
BlkMxL: LOCATE Mx, 1: PRINT Line$: PRINT Line$: LOCATE Mx, 1: RETURN
PauseLine: GOSUB BlkMxL: PRINT "Press any key to continue."
DO: LOOP UNTIL INKEY$ <> Null$: GOSUB BlkMxL
RETURN
GetKey: Digit = 0: DO: C$ = INKEY$: LOOP UNTIL C$ <> Null$
IF ASC(C$) = 13 THEN RETURN
IF ASC(C$) > 47 AND ASC(C$) < 58 THEN Digit = VAL(C$)
Key$ = UCASE$(C$)
RETURN
GetNum2: DOK = 1: GOSUB BlkMxL: PRINT Pmt$; : LINE INPUT a$
IF a$ = Null$ THEN RETURN
FOR Ign = 1 TO LEN(a$)
B$ = MID$(a$, Ign, 1)
IF ASC(B$) < 46 OR ASC(B$) > 57 THEN DOK = 0
IF ASC(B$) = 47 THEN DOK = 0
NEXT Ign
IF DOK = 0 THEN BEEP: GOTO GetNum2
jwh = VAL(a$)
xjwh = VAL(a$) ' For double precision values
RETURN
ResltScr: IF AnaDone$ = NotDone$ THEN RETURN
IF Abort$ <> Null$ THEN LOCATE BotLin, 1: PRINT Abort$
Col = 40: Row = 2
IF Anal = 1 THEN
LOCATE Row, Col: PRINT USING Fmt01$; "f(x) = "; YofX
Row = Row + 1: LOCATE Row, Col: PRINT "for x ="
FOR I = 1 TO Sol: Row = Row + 1: LOCATE Row, Col: PRINT USING Fmt00$; xOFy(I): NEXT I
END IF
IF Anal = 2 THEN
IF (Sol1 + Sol2 + Sol3) > 15 THEN RETURN
LOCATE Row, Col: PRINT "Inflection for x ="
FOR I = 1 TO Sol3: Row = Row + 1: LOCATE Row, Col: PRINT USING Fmt00$; FPx(I): NEXT I
Row = Row + 1: LOCATE Row, Col: PRINT "Local Max, for x ="
FOR I = 1 TO Sol2: Row = Row + 1: LOCATE Row, Col: PRINT USING Fmt00$; LXmx(I): NEXT I
Row = Row + 1: LOCATE Row, Col: PRINT "Local Min, for x ="
FOR I = 1 TO Sol1: Row = Row + 1: LOCATE Row, Col: PRINT USING Fmt00$; LXmn(I): NEXT I
END IF
IF Anal = 3 THEN
LOCATE Row, Col: PRINT "f1(x) = f2(x)"
Row = Row + 1: LOCATE Row, Col: PRINT "for x ="
FOR I = 1 TO Sol: Row = Row + 1: LOCATE Row, Col: PRINT USING Fmt00$; Y1eqY2atX(I): NEXT I
END IF
RETURN
' *****************
' MAIN OPTIONS MENU
' *****************
' 1. Function Use to select function for INVERSE, INFLEC, MAX, MIN
' Toggles through the 5 DEF FUNC routines.
' 2. Domain minimum Define minimum x value, start of iterations
' 3. Domain maximum Define maximum x value, end of iterations
' 4. # of increments Define number of increment between domain end points.
' Incremental value, xi, equals (DoMax - DoMin) / Npts.
' 5. Analysis type Toggles through, INVERSE, MAX-MIN-INFLECTION, INTERSECTION.
' 6. Function value Use to define given range value for INVERSE analysis.
' 7. Tolerance Use to define tolerance for solutions, least significant digit.
' 8. Time limit (secs) Use to define maximum amout of time spent in iteration loop.
' 9. Function 1 Use to select first function for INTERSECTION analysis.
' 0. Function 2 Use to select second function for INTERSECTION analysis.
' A. Do calculations Use to start iteration analysis.
' B. Show progress Toggles between display/no display of program progress.
' C. End program Terminates program.
MO4Scr:
PmtLine$ = "Main Options. "
COLOR Red, 0: PRINT "MAIN OPTIONS:": COLOR 7, 0
PRINT It1$; : PRINT " Function "; : PRINT FDesc0$
PRINT It2$; : PRINT " Domain minimum "; : PRINT USING Fmt00$; DoMin
PRINT It3$; : PRINT " Domain maximum "; : PRINT USING Fmt00$; DoMax
PRINT It4$; : PRINT " # of increments "; : PRINT USING Fmt00$; Npts
PRINT It5$; : PRINT " Analysis type "; : PRINT Anal$(Anal)
PRINT It6$; : PRINT " Function value "; : : PRINT USING Fmt00$; YofX
PRINT It7$; : PRINT " Tolerance "; : : PRINT USING Fmt00$; Tol#
PRINT It8$; : PRINT " Time limit, secs "; : : PRINT AborTime
PRINT It9$; : PRINT " Function 1 "; : : PRINT FDesc1$
PRINT It0$; : PRINT " Function 2 "; : : PRINT FDesc2$
PRINT ItA$; : PRINT " Do calculations. "; AnaDone$
PRINT ItB$; : PRINT " Progress "; ProgStat$
PRINT ItC$; : PRINT " End program "
RETURN
' ************************
' END OF MAIN OPTIONS MENU
' ************************
Options:
MO4: CLS : GOSUB MO4Scr: GOSUB ResltScr
GOSUB BlkMxL
PRINT "Choose from "; : PRINT PmtLine$
GOSUB GetKey: MO4Key$ = Key$
IF Key$ <> ItF$ AND Key$ <> ItG$ AND Key$ <> ItH$ THEN AnaDone$ = NotDone$
SELECT CASE MO4Key$
CASE It1$
Feq0 = Feq0 + 1: IF Feq0 > FeqMax THEN Feq0 = 1
x = .0001235: Feq = Feq0: GOSUB Selfunc
FDesc0$ = FDesc$
CASE It2$
Pmt$ = "Enter minimum x value of domain. ": jwh = DoMin
GOSUB GetNum2: DoMin = jwh
CASE It3$
Pmt$ = "Enter maximum x value of domain. ": jwh = DoMax
GOSUB GetNum2: DoMax = jwh
CASE It4$
Pmt$ = "Enter number of domain increments. ": jwh = Npts
GOSUB GetNum2: Npts = jwh
CASE It5$
Anal = Anal + 1: IF Anal > 3 THEN Anal = 1
CASE It6$
Pmt$ = "Enter function value. ": jwh = YofX
GOSUB GetNum2: YofX = jwh
CASE It7$
Pmt$ = "Enter the tolerance for solution. ": xjwh = Tol#
GOSUB GetNum2: Tol# = xjwh
CASE It8$
Pmt$ = "Enter abort analysis time limit, secs. ": jwh = AborTime
GOSUB GetNum2: AborTime = jwh
CASE It9$
Feq1 = Feq1 + 1: IF Feq1 > 5 THEN Feq1 = 1
x = .0000123: Feq = Feq1: GOSUB Selfunc: FDesc1$ = FDesc$
CASE It0$
Feq2 = Feq2 + 1: IF Feq2 > 5 THEN Feq2 = 1
x = .0000123: Feq = Feq2: GOSUB Selfunc: FDesc2$ = FDesc$
CASE ItA$
GOSUB Calculations
CASE ItB$
Progress = -Progress
IF Progress = 1 THEN ProgStat$ = "On" ELSE ProgStat$ = "Off"
CASE ItC$
CLS : STOP
END SELECT
GOTO MO4
' End of MAIN PROGRAM: Data entry and change, results display.
Selfunc: ' Select the required function according to Feq and evaluate f(x).
' Define FUNCTIONs using QBASIC software development.
IF Feq = 1 THEN y = FUNC1#(x, FDesc$)
IF Feq = 2 THEN y = FUNC2#(x, FDesc$)
IF Feq = 3 THEN y = FUNC3#(x, FDesc$)
IF Feq = 4 THEN y = FUNC4#(x, FDesc$)
IF Feq = 5 THEN y = FUNC5#(x, FDesc$)
RETURN
' Iteration procedures follow.
Calculations:
GOSUB BlkMxL: PRINT "Press any key to abort."
' Use number of domain intervals to find initial x increment value.
Pdx = (DoMax - DoMin) / Npts
Analysis:
GOSUB BlkMxL: PRINT "Press any key to abort."
SELECT CASE Anal
CASE 1 ' Inverse
Feq = Feq0 ' Set select function variable
Sol = 0 ' Reset solution counter
x = DoMin - Pdx ' Incremention at top of loop, back up one step
xi = Pdx ' Set increment to initial value
TimeStop = TIMER + AborTime ' Set abort timer
Abort$ = Null$ ' Set abort message to null string
DO:
x = x + xi ' Increment x
GOSUB Selfunc ' Evaluate selected function at x
y1 = y ' Retain y at this point
x = x + xi: GOSUB Selfunc ' Repeat for next adjacent x value
y2 = y: x = x - xi
' For positive slope, test if given value, YofX, is between these two values,
' if it is then back up one step, refine increment by .1*xi, back up one step
' of refined increment.
IF YofX >= y1 AND YofX < y2 THEN GOSUB BackUp
' For negative slope, repeat same process of backing up.
IF YofX <= y1 AND YofX > y2 THEN GOSUB BackUp
IF xi < (Tol# * .001) THEN ' Assume discontinuity,
x = x + Pdx: xi = Pdx ' Skip this interval and reset increment value
GOTO SkipPoint1 ' Skip over rest of loop.
END IF
IF ABS(YofX - y1) <= Tol# THEN ' If absolute difference is within tolerance
Sol = Sol + 1 ' increment solution counter
xOFy(Sol) = x ' save x solution
xi = Pdx ' reset increment to original value
x = x + xi ' step forward one to avoid repeating
END IF ' this solution
IF Progress = 1 THEN ' If progress switch is on then display this
FOR I = 1 TO 10: LOCATE I, 40: PRINT BlkHL$: NEXT I
LOCATE 1, 40: PRINT " Progress"
LOCATE 2, 40: PRINT "xinc = "; xi
LOCATE 3, 40: PRINT " x = "; x
LOCATE 4, 40: PRINT " y1 = "; y1
LOCATE 5, 40: PRINT "YofX = "; YofX
LOCATE 6, 40: PRINT " y2 = "; y2
LOCATE 7, 40: PRINT " Sol = "; Sol
PauseTime = TIMER + .05 ' Hesitate display to avoid flicking
DO: LOOP UNTIL TIMER > PauseTime
END IF
SkipPoint1: ' Jump here if function discontinuity
GOSUB ChkTime ' Check to see if time limit has been exceeded,
' if true the set Abort$ message to display such.
' Do loop until one of these has occurred,
' INKEY$ is manual abort
LOOP UNTIL x >= DoMax OR Sol >= 10 OR TIMER > TimeStop OR INKEY$ <> ""
CASE 2 ' Inflextion, local max, local min points
Feq = Feq0
xi = Pdx
TimeStop = TIMER + AborTime
Abort$ = Null$
Slope = 1 ' Slope switch: 1 = positive, -1 = negative
Rate = 1 ' Slope rate: 1 = increasing, -1 = decreasing
x = DoMin ' Evaluate f(x) at initial point and next 2 points
GOSUB Selfunc: y1 = y
x = x + xi
GOSUB Selfunc: y2 = y
x = x + xi
GOSUB Selfunc: y3 = y
dydx1 = (y2 - y1) / xi ' Slope over first interval
dydx2 = (y3 - y2) / xi ' Slope over second interval
IF dydx1 < 0 THEN Slope = -1 ' Reset slope switch if negative
IF dydx1 > dydx2 THEN Rate = -1 ' Reset rate switch if decreasing
x = DoMin - xi ' Incremention at top of loop
Sol1 = 0: Sol2 = 0: Sol3 = 0 ' Reset solution counters
LocMin = -1: LocMax = -1 ' Set local minimum, maximum switches
Inflex1 = -1: Inflex2 = -1 ' Set inflection switches
DO
x1 = x ' Retain current point as a solution
GOSUB Selfunc: y1 = y ' Evaluate f(x) at current point and
x = x + xi: GOSUB Selfunc: y2 = y ' the next two points
x = x + xi: GOSUB Selfunc: y3 = y
x = x - xi ' Back up one point for next iteration
dydx1# = (y2 - y1) / xi ' Slope of first interval
dydx2# = (y3 - y2) / xi ' Slope of second interval
' If current slope is negative and the previous slope positive then
' a local maximum point present. Set LocMax switch and do the BackUp.
' If current slope is positive and the previous slope negative then
' a local minimum point is present. Set LocMin switch and do the BackUp.
IF dydx1# < 0 AND Slope = 1 THEN LocMax = 1: GOSUB BackUp ' + to - change
IF dydx1# > 0 AND Slope = -1 THEN LocMin = 1: GOSUB BackUp ' - to + change
' If current slope rate is decreasing and rate was increasing or
' if current slope rate is increasing and rate was decreasing then
' an inflection point is present. Set Inflec1 switch and do BackUp routine.
IF dydx1# > dydx2# AND Rate = 1 THEN Inflex1 = 1: GOSUB BackUp ' - to +
IF dydx1# < dydx2# AND Rate = -1 THEN Inflex2 = 1: GOSUB BackUp ' + to -
IF xi < Tol# THEN ' If accuracy has been achieved then
IF LocMin = 1 AND Sol1 < 10 THEN ' determine which type of solution it is
Sol1 = Sol1 + 1 ' increment the solution counter
LYmn(Sol1) = y1: LXmn(Sol1) = x1 ' save the y and x values
LocMin = -1: Slope = 1 ' reset switches appropriately
END IF ' also check for maximum number
IF LocMax = 1 AND Sol2 < 10 THEN ' of solutions in arrays
Sol2 = Sol2 + 1
LYmx(Sol2) = y1: LXmx(Sol2) = x1
LocMax = -1: Slope = -1
END IF
IF (Inflex1 = 1 OR Inflex2 = 1) AND Sol3 < 10 THEN
Sol3 = Sol3 + 1: FPy(Sol3) = y1: FPx(Sol3) = x1
IF Inflex1 = 1 THEN Rate = -1
IF Inflex2 = 1 THEN Rate = 1
Inflex1 = -1: Inflex2 = -1
END IF
xi = Pdx
END IF
IF Progress = 1 THEN ' If Progress switch is one then display this
FOR I = 1 TO 10: LOCATE I, 40: PRINT BlkHL$: NEXT I
LOCATE 1, 40: PRINT " STATUS"
LOCATE 2, 40: PRINT "xinc = "; xi
LOCATE 3, 40: PRINT " x = "; x
LOCATE 4, 40: PRINT " y1 = "; y1
LOCATE 5, 40: PRINT "Solutions:"
LOCATE 6, 40: PRINT "LocMin = "; Sol1
LOCATE 7, 40: PRINT "LocMax = "; Sol2
LOCATE 8, 40: PRINT "Inflec = "; Sol3
PauseTime = TIMER + .05 ' Hesitate to avoid flicking
DO: LOOP UNTIL TIMER > PauseTime
END IF
GOSUB ChkTime
LOOP UNTIL x > DoMax OR TIMER > TimeStop OR INKEY$ <> Null$
CASE 3 ' Intersection
Sol = 0 ' Reset solution counter
x = DoMin - Pdx ' Incremention at top of loop
xi = Pdx ' Set initial increment value
Test = -1 ' Test switch: 1st function > 2nd function = 1
' 1st function < 2nd function = -1
TimeStop = TIMER + AborTime ' Set loop timeout limit
Abort$ = Null$ ' and null abort message
Feq = Feq1: GOSUB Selfunc: y1 = y ' Evaluate 1st function
Feq = Feq2: GOSUB Selfunc: y2 = y ' Evaluate 2nd function
IF y1 > y2 THEN Test = 1 ' Reset Test switch if required
DO
x = x + xi ' Increment to next x value
Feq = Feq1: GOSUB Selfunc: y1 = y ' Evaluate 1st function
Feq = Feq2: GOSUB Selfunc: y2 = y ' Evaluate 2nd function
' If increment value is now less than tolerance
' assume a discontinuity at this interval and skip over it.
IF xi < (Tol# * .001) THEN
Test = -Test ' For discontinuity relative position has changed
x = x + Pdx: xi = Pdx ' Skip this interval and reset increment value
GOTO SkipPoint3 ' Skip over rest of loop.
END IF
' If relative position of functions has changed then do BackUp routine
IF y1 > y2 AND Test = -1 THEN GOSUB BackUp ' 1st > 2nd and was 1st < 2nd
IF y1 < y2 AND Test = 1 THEN GOSUB BackUp ' 1st < 2nd and was 1st > 2nd
IF ABS(y1 - y2) < Tol# THEN ' If accuracy has been achieved then
Sol = Sol + 1 ' increment solution counter
Y1eqY2atX(Sol) = x ' save x solution and
y1atx(Sol) = y1 ' use 1st function value as y solution
Test = -Test ' reset Test switch
xi = Pdx ' reset increment to original value
END IF
IF Progress = 1 THEN ' If Progress switch is on then display this
FOR I = 1 TO 10: LOCATE I, 40: PRINT BlkHL$: NEXT I
LOCATE 1, 40: PRINT " STATUS"
LOCATE 2, 40: PRINT "xinc = "; xi
LOCATE 3, 40: PRINT " x = "; x
LOCATE 4, 40: PRINT "f1(x)= "; y1
LOCATE 5, 40: PRINT "f2(x)= "; y2
LOCATE 6, 40: PRINT " Sol = "; Sol
PauseTime = TIMER + .05 ' Hesitate to avoid flickering
DO: LOOP UNTIL TIMER > PauseTime
END IF
SkipPoint3: ' Jump here if discontiuity this inteveral.
GOSUB ChkTime
LOOP UNTIL x > DoMax OR Sol >= 10 OR TIMER > TimeStop OR INKEY$ <> Null$
END SELECT
AnaDone$ = Done$
RETURN
BackUp:
IF xi > (Tol# * .001#) THEN x = x - 1 * xi: xi = xi * .1: x = x - xi
RETURN
ChkTime:
IF TIMER > TimeStop THEN
BEEP: Abort$ = "Exceeded time limit, analysis incomplete."
END IF
RETURN
' ********************
' FUNCTION DEFINITIONS
' ********************
'
' User should use QB45 or other basic language environment to define the
' FUNCTIONS FUNC1#(x, D$), FUNC2#(x, D$), ... which are to be analyized.
' The code can be simple equations or complex numerical models within each
' function definition.
'
' There are probably several ways to manage these function definitions for
' different applications or problems. Two suggested methods follow.
'
' First method.
' 1. Save the MAIN PROGRAM code, ie, everything but FUNCTIONS under one file
' name. SAVE AS PLAIN TEXT ASCII.
' 2. Develop function definitions and save under a filename unique for each
' problem. SAVE AS PLAIN TEXT ASCII.
' 3. Append the FUNCTIONs file to the MAIN PROGRAM file and run the program.
' Second method.
' 1. Save MAIN PROGRAM and FUNCTIONs definitions as one file for each
' application or problem.
' 2. Use of this methods allows easier development of FUNCTION code within
' the language environment.
END ' End of MAIN PROGRAM code
FUNCTION FUNC1# (x, D$)
D$ = "sin(x)"
FUNC1# = 1! * SIN(x)
END FUNCTION
FUNCTION FUNC2# (x, D$)
D$ = "tan(x)"
FUNC2# = TAN(x)
END FUNCTION
FUNCTION FUNC3# (x, D$)
D$ = "y = .75"
FUNC3# = .75
END FUNCTION
FUNCTION FUNC4# (x, D$)
D$ = "x^2+x+0"
FUNC4# = x ^ 2 + x + 0
END FUNCTION
FUNCTION FUNC5# (x, D$)
D$ = "y=.5"
FUNC5# = .5
END FUNCTION
FUNCTION LGT (x#)
LGT = LOG(x#) / LOG(10#)
END FUNCTION
' Delete the below lines and this one.