home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / MBUG / MBUG058.ARC / CURVFIT.BAS < prev    next >
BASIC Source File  |  1979-12-31  |  7KB  |  243 lines
  1. 10 REM
  2. 20 REM   POLYNOMIAL  CURVFIT
  3. 30 REM   By Ralph Johnson 1-1982
  4. 35 REM
  5. 40       REM  MAIN PROGRAM
  6. 100 GOSUB 1000     REM SET UP ARRAYS
  7. 110 GOSUB 1500  REM INPUT DATA
  8. 120 GOSUB 2000  REM CLEAR AND LOAD ARRAYS
  9. 130 GOSUB 2500  REM SIMUL. EQS. SOLUTION
  10. 135 GOSUB 3000  REM SUM OF ERRORS SQRD
  11. 140 GOSUB 3500  REM OUTPUT
  12. 150 GOTO 4000   REM DECISIONS
  13. 160 GOSUB 4500  REM EVALUATE
  14. 170 GOSUB 5000  REM PLOT
  15. 180 GOTO 4000
  16. 190    REM END MAIN
  17. 200 REM
  18. 210    REM SET UP ARRAY MEMORY
  19. 220 REM
  20. 1000 DIM A(20,20),B(20),C(20),X(20),Y(20)
  21. 1010 DIM XX(120),YY(120),Q(45),YZ(20)
  22. 1020 REM    GREETING
  23. 1025 PRINT CHR$(&H1A)
  24. 1030 PRINT "******* POLYNOMIAL  CURVE  FIT *******"
  25. 1040 PRINT:PRINT
  26. 1050 PRINT"    THIS PROGRAM PERFORMS A LEAST SQUARED"
  27. 1060 PRINT"CURVE FIT FOR A GENERAL POLYNOMIAL IN X AND Y"
  28. 1070 PRINT"IT WILL ALSO EVALUATE AND PLOT THE FUNCTION"
  29. 1080 PRINT"IF YOU WISH"
  30. 1090 PRINT:PRINT:PRINT
  31. 1100 RETURN
  32. 1110    REM END SETUP MEMORY
  33. 1500    REM INPUT SUBROUTINE 
  34. 1510    REM RESET ARRAYS
  35. 1520 FOR I%=0 TO 20 
  36. 1530    B(I%)=0
  37. 1540    Q(I%)=0
  38. 1560 NEXT I%
  39. 1565 IF Z% =2 THEN GOTO 1670
  40. 1570 INPUT"HOW MANY DATA PTS. ARE THERE";DATAPTS%
  41. 1575 PRINT CHR$(&H1A)
  42. 1580 FOR I%=1 TO 5
  43. 1590    PRINT
  44. 1600 NEXT I%
  45. 1610 PRINT"NOW INPUT THE DATA.":PRINT:PRINT
  46. 1620 FOR I%=1 TO DATAPTS%
  47. 1630    PRINT"X(";I%;"), Y(";I%;")";
  48. 1640    INPUT X(I%),Y(I%)
  49. 1650 NEXT I%
  50. 1660 REM
  51. 1670 INPUT"WHAT ORDER FIT WOULD YOU LIKE";ORDER%
  52. 1680     N%=ORDER%+1
  53. 1690 REM 
  54. 1700    REM VALIDITY CHECK
  55. 1710 REM
  56. 1720 IF ORDER%>=DATAPTS% THEN PRINT"NOT ENOUGH DATA FOR THAT ORDER FIT": GOTO 1670
  57. 1730 PRINT CHR$(&H1A)
  58. 1750 RETURN
  59. 1760    REM END INPUT
  60. 1970    REM MATRIX LOAD
  61. 1980  REM GENERATE VALUES FOR ARRAY
  62. 2000 FOR I%=0 TO ORDER%*2
  63. 2010    FOR J%=1 TO DATAPTS%
  64. 2020        Q(I%)=X(J%)^I%+Q(I%)
  65. 2030    NEXT J%
  66. 2040    Q(I%)=2*Q(I%)
  67. 2050 NEXT I%
  68. 2060 REM LOAD
  69. 2070 FOR I%=0 TO ORDER%
  70. 2080    FOR J%=0 TO ORDER%
  71. 2090        A(I%+1,J%+1)=Q(I%+J%)
  72. 2100    NEXT J%
  73. 2110 NEXT I%
  74. 2120 FOR I%=0 TO ORDER%
  75. 2130    FOR J%=1 TO DATAPTS%
  76. 2140        B(I%+1)=2*X(J%)^I%*Y(J%)+B(I%+1)
  77. 2150    NEXT J%
  78. 2160 NEXT I%
  79. 2170 RETURN
  80. 2180    REM END MATRIX LOAD
  81. 2500    REM SIMULTANEOUS EQU SOLN
  82. 2510 FLAG%=0
  83. 2520  IF N%<>0 THEN GOTO 2560
  84. 2530  IF A(1,1)=0 THEN FLAG%=1: RETURN
  85. 2540  X(1)=B(1)/A(1,1)
  86. 2550  RETURN
  87. 2560  M%=N%-1
  88. 2570  FOR I%=1 TO M%
  89. 2580    H=ABS(A(I%,I%))
  90. 2590    L%=I%
  91. 2600    I1%=I%+1
  92. 2610    FOR J%=I1% TO N%
  93. 2620       IF ABS(A(J%,I%)) < H THEN 2650
  94. 2630       H=ABS(A(J%,I%))
  95. 2640       L%=J%
  96. 2650    NEXT J%
  97. 2660    IF H=0 THEN FLAG%=1: RETURN
  98. 2670    IF L%=I% THEN 2760
  99. 2680    FOR J%=1 TO N%
  100. 2690       G=A(L%,J%)
  101. 2700       A(L%,J%)=A(I%,J%)
  102. 2710       A(I%,J%)=G
  103. 2720    NEXT J%
  104. 2730    G=B(L%)
  105. 2740    B(L%)=B(I%)
  106. 2750    B(I%)=G
  107. 2760    FOR J%=I1% TO N%
  108. 2770       T=A(J%,I%)/A(I%,I%)
  109. 2780       FOR K%=I1% TO N%
  110. 2790        A(J%,K%)=A(J%,K%)-T*A(I%,K%)
  111. 2800       NEXT K%
  112. 2810       B(J%)=B(J%)-T*B(I%)
  113. 2820    NEXT J%
  114. 2830  NEXT I%
  115. 2840 IF A(N%,N%)=0 THEN FLAG%=1: RETURN
  116. 2850 C(N%)=B(N%)/A(N%,N%)
  117. 2860 I%=N%-1
  118. 2870 S=0
  119. 2880 I1%=I%+1
  120. 2890 FOR J%=I1% TO N%
  121. 2900    S=S+A(I%,J%)*C(J%)
  122. 2910 NEXT J%
  123. 2920 C(I%)=(B(I%)-S)/A(I%,I%)
  124. 2930 I%=I%-1
  125. 2940 IF I%>0 THEN 2870
  126. 2950 RETURN
  127. 2960     REM END SIMULTANEOUS EQU
  128. 3000     REM SUM OF ERRORS SQRD
  129. 3010 FOR I%=1 TO DATAPTS%
  130. 3020     YZ(I%)=0
  131. 3030 NEXT I%
  132. 3040  E=0
  133. 3050  FOR I%=1 TO DATAPTS%
  134. 3060     FOR J%=1 TO N%
  135. 3070        YZ(I%)=C(J%)*X(I%)^(J%-1)+YZ(I%)
  136. 3080     NEXT J%
  137. 3090     E=E+(Y(I%)-YZ(I%))^2
  138. 3100 NEXT I%
  139. 3110 RETURN
  140. 3120 REM END SUM ERRORS SQRD
  141. 3500     REM OUTPUT STATEMENTS
  142. 3510 IF FLAG%=1 THEN PRINT "AMBIGUOUS DATA. NO SOLUTION FOUND":RETURN
  143. 3520 PRINT"LIST OF COEFFICIENTS":PRINT
  144. 3530 PRINT"POWER OF X","COEFFICIENT":PRINT
  145. 3540 FOR I%=N% TO 1 STEP -1
  146. 3550     PRINT "X^";I%-1,C(I%)
  147. 3560 NEXT I%
  148. 3570 PRINT "SUM OF ERRORS SQRD=";E
  149. 3580 PRINT
  150. 3590 FOR I%=1 TO DATAPTS%
  151. 3600    PRINT"X=";X(I%),"Y=";YZ(I%),"DESIRED Y:";Y(I%)
  152. 3610 NEXT I%
  153. 3620    REM Printer
  154. 3630 PRINT:PRINT
  155. 3700 INPUT "Do you want a print out of this (y/n)";PR$
  156. 3710 IF PR$<>"y" AND PR$<>"Y" THEN GOTO 3830
  157. 3720 LPRINT "LIST OF COEFFICIENTS": LPRINT " "
  158. 3730 LPRINT " ": LPRINT "POWER OF X", "COEFFICIENT": LPRINT" "
  159. 3740 FOR I%=N% TO 1 STEP -1
  160. 3750    LPRINT "X^";I%-1,C(I%)
  161. 3760 NEXT I%
  162. 3770 LPRINT : LPRINT "SUM OF ERRORS SQRD=";E
  163. 3780 LPRINT" "
  164. 3790 FOR I%=1 TO DATAPTS%
  165. 3800    LPRINT "x=";X(I%),"y=";YZ(I%),"Desired y:";Y(I%)
  166. 3810 NEXT I%
  167. 3820 RETURN
  168. 3830 REM end output
  169. 4000 REM OPTIONS
  170. 4010 REM
  171. 4020 INPUT"HIT ENTER TO CONTINUE";A$
  172. 4030 PRINT CHR$(&H1A): PRINT: PRINT
  173. 4200 PRINT"Enter the desired option": PRINT: PRINT
  174. 4210 PRINT"    1)EVALUATE THE FUNCTION"
  175. 4220 PRINT"    2)TRY A DIFFERENT ORDER FIT"
  176. 4230 PRINT"    3)TRY A NEW DATA SET"
  177. 4235 PRINT"    4)EXIT THE PROGRAM"
  178. 4240 PRINT:INPUT"Enter your choice";Z%
  179. 4250 ON Z% GOTO 160,110,110,4260
  180. 4260 STOP
  181. 4270 REM    END OPTIONS
  182. 4500      REM EVALUATE
  183. 4520 INPUT"SELECT THE RANGE YOU WANT TO EVALUATE min,max";XINIT,XFIN
  184. 4530 INCR=(XFIN-XINIT)/120
  185. 4540 ZO=0
  186. 4550 FOR I%=0 TO 120
  187. 4560     XX(I%)=XINIT+INCR*I%
  188. 4570     YY(I%)=0
  189. 4580    FOR J%=1 TO N%
  190. 4590        YY(I%)=C(J%)*XX(I%)^(J%-1)+YY(I%)
  191. 4600    NEXT J%
  192. 4610    IF I%=0 THEN YMIN=YY(0):YMAX=YMIN
  193. 4620    IF YY(I%)>YMAX THEN YMAX=YY(I%)
  194. 4630    IF YY(I%)<YMIN THEN YMIN=YY(I%)
  195. 4640 NEXT I%
  196. 4650    REM output evaluate
  197. 4680 PRINT CHR$(&H1A)
  198. 4690 PRINT"X AND Y VALUES":PRINT
  199. 4700 PRINT"X","Y","X","Y"
  200. 4710 FOR I%=0 TO 60 STEP 6
  201. 4720    PRINT XX(I%),YY(I%),XX(I%+60),YY(I%+60)
  202. 4730 NEXT I%
  203. 4740 PRINT
  204. 4750 PRINT"Y MAX.:";YMAX,"Y MIN.:";YMIN
  205. 4760 INPUT "DO YOU WANT THIS PRINTED OUT (y/n)";PR$
  206. 4770 IF PR$<>"Y" AND PR$<>"y" THEN RETURN
  207. 4790 LPRINT"X and Y values" : LPRINT " "
  208. 4800 LPRINT"X","Y","X","Y"
  209. 4810 FOR I%=0 TO 60 STEP 6
  210. 4820    LPRINT XX(I%),YY(I%),XX(I%+60),YY(I%+60)
  211. 4830 NEXT I%
  212. 4840 LPRINT" "
  213. 4850 LPRINT "Y MIN=";YMIN,"Y MAX=";YMAX
  214. 4860 RETURN
  215. 5000    REM PLOT
  216. 5010 INPUT"DO YOU WANT A PLOT OF THE FUNCTION (y/n)";PLOT$
  217. 5020 IF PLOT$="N" OR PLOT$="n" THEN RETURN
  218. 5030 INPUT"SELECT THE RANGE OF Y (ymin,ymax)";
  219.       YMIN,YMAX
  220. 5040 INPUT"SELECT THE COLUMN WIDTH";WDTH 
  221. 5045 YRANGE=YMAX-YMIN
  222. 5050 ZERO=(0-YMIN)*WDTH/YRANGE+.5
  223. 5060 YRANGE=YMAX-YMIN
  224. 5070 FOR I%=0 TO 119 STEP 2
  225. 5080    PLOT$=" "
  226. 5090    YPLOT1%=INT((YY(I%)-YMIN)*WDTH/YRANGE+.5)
  227. 5100    YPLOT2%=INT((YY(I%+1)-YMIN)*WDTH/YRANGE+.5)
  228. 5110    FOR J%= 1 TO WDTH
  229. 5120       IF YPLOT1%=J% AND YPLOT2%=J% THEN
  230.         PLOT$=PLOT$+";": GOTO 5160
  231. 5130       IF YPLOT1%=J% THEN PLOT$=PLOT$+"'":
  232.              GOTO 5160
  233. 5140       IF YPLOT2%=J% THEN PLOT$=PLOT$+",":
  234.          GOTO 5160
  235. 5150       IF INT(ZERO)=J% THEN PLOT$=PLOT$+"|":
  236.         GOTO 5160
  237. 5155       PLOT$=PLOT$+" "
  238. 5160    NEXT J%
  239. 5170    PRINT PLOT$
  240. 5175    LPRINT PLOT$
  241. 5180 NEXT I%
  242. 5190 RETURN
  243.