home *** CD-ROM | disk | FTP | other *** search
/ 64'er 1987 December / 64er_Magazin_87-12_1987_Markt__Technik_de_Side_A.d64 / polynome.c64 (.txt) < prev    next >
Encoding:
Commodore BASIC  |  1987-01-01  |  7.6 KB  |  226 lines
  1. 1 rem ----- nachladen grafik -------
  2. 2 if a$="n" then 100
  3. 3 printchr$(147)chr$(17)"grafik nachladen (j/n)"
  4. 4 get a$:if a$="" then 4
  5. 5 if a$="n" then 100
  6. 6 a$="n":load"sysgraf.obj",8,1
  7. 10 rem*********************************
  8. 20 rem*                               *
  9. 30 rem*       kurvenanpassung         *
  10. 40 rem*                               *
  11. 50 rem*    statistik-programm zur     *
  12. 60 rem*  regressionsanalyse mittels   *
  13. 65 rem*     polynomen incl.grafik     *
  14. 70 rem*  heimo ponnath  hamburg 1987  *
  15. 80 rem*       c64  - version          *
  16. 90 rem*********************************
  17. 100 clr:sys 49152:sys 49242:rem grafikspeicher sichern
  18. 105 rem graphic1,1:gosub 4000
  19. 110 rem ----- variable ----------------
  20. 120 m=0:s=0:n=0:i=0:j=0:k=0:g=0:a=0:g2=0:zz=0:hh=0:q=0:p=0:a0=0:b=0
  21. 130 bb=319:bh=199:w=0:mx=-1e12:lx=1e12:my=-1e12:ly=1e12:r=0
  22. 140 xu=0:xo=0:yu=0:yo=0:ra=0:rd=0:ta=0:tb=0:x=0:y=0:x1=0:y1=0
  23. 150 a$="":b$=""
  24. 160 rem
  25. 170 deffn g(n)=int(.2*sqr(5*(a-10*n-17))-4):rem max. polynomgrad berechnen
  26. 180 deffn x(x)=ra*x+ta:deffn y(y)=rd*y+tb
  27. 190 rem
  28. 200 rem ----- titel,erklaerung --------
  29. 210 poke 53280,0:poke 53281,0:printchr$(30)
  30. 220 printchr$(147)chr$(18)"            polynomanpassung            "chr$(146)
  31. 230 print
  32. 240 print" durch eine anzahl von n punkten aus"
  33. 250 print"wertepaaren legt dieses programm das am"
  34. 260 print"besten angepasste polynom der form"
  35. 270 print"y=a0+a1*x+a2*x^2+a3*x^3+... . der grad"
  36. 275 print"des polynoms ist frei waehlbar. der"
  37. 277 print"korrelationskoeffizient r und die"
  38. 280 print"standardabweichung s werden angegeben"
  39. 290 print"und sie koennen beliebige y-werte"
  40. 300 print"aus eingegebenen x-werten berechnen.":print
  41. 310 print" ein scatterdiagramm und die ermittelte"
  42. 320 print"kurve werden gezeichnet. auf diese"
  43. 330 print"weise kann die qualitaet der anpassung"
  44. 340 print"eingeschaetzt werden.":print
  45. 350 print" sogenannte ausreisser-werte sollten"
  46. 360 print"vor einer genaueren berechnung noch"
  47. 370 print"entfernt werden.":print:print
  48. 380 printchr$(18)"taste druecken!"chr$(146)
  49. 390 get a$:if a$="" then 390
  50. 400 rem ----- hauptmenue --------------
  51. 410 printchr$(147):print:print:print:print
  52. 420 printtab(4)"werte von hand eingeben.......1":print
  53. 430 printtab(4)"werte aus datei lesen.........2":print
  54. 440 printtab(4)"grafik zeigen.................3":print
  55. 450 printtab(4)"textmodus einschalten.........4":print
  56. 460 printtab(4)"polynomfunktion berechnen.....5":print
  57. 470 printtab(4)"werte berechnen...............6":print
  58. 480 printtab(4)"programmende..................7":print:print
  59. 490 printtab(10)chr$(18)"bitte waehlen sie!"chr$(146)
  60. 500 get a$:if val(a$)<1 or val(a$)>7 then 500
  61. 510 printchr$(147):if val(a$)=7 then end
  62. 520 on val(a$) gosub 1000,2000,3000,4000,5000,6000
  63. 530 goto 410
  64. 540 rem ----- ende hauptprogramm ------
  65. 1000 rem ----- werte von hand ---------
  66. 1005 gosub 4000:rem textmodus
  67. 1010 if w=1 then print"werte schon vorhanden!":for i=0 to 500:next i:return
  68. 1020 w=1
  69. 1030 print"wieviele werte werden verwendet ?":inputn:print
  70. 1040 dim w(1,n)
  71. 1043 a=fre(0)-2000:rem freier speicherplatz c64
  72. 1045 rem a=fre(1)-2000
  73. 1047 g=fng(n):rem maximaler polynomgrad
  74. 1050 print"bitte wertepaare eingeben!":print
  75. 1060 for i=1 to n
  76. 1070 printi,"x=";:inputw(0,i):printchr$(145),,"y=";:inputw(1,i):print
  77. 1080 gosub 1300:rem zwischenwerte berechnen
  78. 1090 next i
  79. 1100 printchr$(147):print:print"sollen die werte gespeichert werden?"
  80. 1110 get a$:if a$<>"j" and a$<>"n" then 1110
  81. 1120 if a$="n" then 1190
  82. 1130 print:print"name der datei (11 zeichen)";:input b$
  83. 1140 b$=left$(b$,11)+".dat"+",s,w"
  84. 1150 open 1,8,2,b$
  85. 1160 print#1,n
  86. 1170 for i=1 to n:print#1,w(0,i):print#1,w(1,i):next i
  87. 1180 close 1
  88. 1190 gosub 1500:rem scatterdiagramm zeichnen
  89. 1200 return
  90. 1300 rem --- zwischenwerte berechnen --
  91. 1360 if w(0,i)>mx then mx=w(0,i)
  92. 1370 if w(0,i)<lx then lx=w(0,i)
  93. 1380 if w(1,i)>my then my=w(1,i)
  94. 1390 if w(1,i)<ly then ly=w(1,i)
  95. 1400 return
  96. 1500 rem --- scatterdiagramm ----------
  97. 1505 dim a(2*g+1),r(g+1,g+2),t(g+2):rem arrays fuer berechnungen
  98. 1510 for i=1 to n-1:rem sortieren nach x
  99. 1520 for j=i+1 to n
  100. 1530 if w(0,i)<=w(0,j) then 1560
  101. 1540 w(0,0)=w(0,i):w(1,0)=w(1,i):w(0,i)=w(0,j):w(1,i)=w(1,j)
  102. 1550 w(0,j)=w(0,0):w(1,j)=w(1,0)
  103. 1560 next j:next i
  104. 1570 sys 49152:sys 49180:sys 49202,6,0:rem grafik loeschen farbe
  105. 1580 rem graphic1,1:color0,1:color1,7
  106. 1590 sys 49352,0,0,319,0,1:sys 49352,319,0,319,199,1
  107. 1600 rem draw 1,0,0 to 319,0 to 319,199 to 0,199 to 0,0
  108. 1610 sys 49352,319,199,0,199,1:sys 49352,0,199,0,0,1:rem rahmen
  109. 1620 xu=lx-(mx-lx)*.02:xo=mx+(mx-lx)*.02
  110. 1630 yu=ly-(my-ly)*.02:yo=my+(my-ly)*.02
  111. 1640 ra=bb/(xo-xu):rd=-bh/(yo-yu)
  112. 1650 ta=-bb*xu/(xo-xu):tb=bh*yo/(yo-yu)
  113. 1660 for i=1 to n
  114. 1670 x=fnx(w(0,i)):y=fny(w(1,i))
  115. 1680 sys49352,x-3,y,x+3,y,1:sys49352,x,y-3,x,y+3,1:rem kreuz
  116. 1681 rem draw1,x-3,y to x+3,y:draw1,x,y-3 to x,y+3
  117. 1690 next i
  118. 1700 get a$:if a$="" then 1700
  119. 1710 sys 49242:rem textmodus
  120. 1711 rem if peek(238)=79 then graphic5:else graphic0
  121. 1720 print"xu =  "lx,"xo =  "mx"
  122. 1730 [153]"yu =  "ly,"yo =  "my"
  123. 1740 get a$:if a$="" then 1740
  124. 1750 return
  125. 2000 rem ----- werte aus datei --------
  126. 2005 gosub 4000:rem textmodus
  127. 2010 if w=1 then print"werte schon vorhanden!":for i=0 to 500:next i:return
  128. 2020 w=1
  129. 2030 print" die datei muss ein bestimmtes format"
  130. 2040 print"haben:     1.anzahl der wertepaare"
  131. 2050 print"           1.wert x, 1.wert y"
  132. 2060 print"           2.wert x, 2.wert y ...":print
  133. 2070 print"diese dateien werden unter menuepunkt 1"
  134. 2080 print"erstellt. sie tragen die endung .dat .":print
  135. 2090 print" alles klar..1  ach soo..2"
  136. 2100 get a$:if val(a$)<1 or val(a$)>2 then 2100
  137. 2110 if val(a$)=2 then w=0:return
  138. 2120 print:print"wie heisst denn die datei (endung .dat)"
  139. 2130 input b$
  140. 2140 b$=b$+",s,r"
  141. 2150 open1,8,2,b$
  142. 2160 input#1,n
  143. 2170 dim w(1,n)
  144. 2173 a=fre(0)-2000:rem freier speicherplatz c64
  145. 2175 rem a=fre(1)-2000
  146. 2177 g=fng(n):rem maximaler polynomgrad
  147. 2180 for i=1 to n
  148. 2190 input#1,w(0,i):input#1,w(1,i)
  149. 2200 gosub 1300:rem zwischenwerte berechnen
  150. 2210 next i
  151. 2220 close 1
  152. 2230 gosub 1500:rem scatterdiagramm
  153. 2240 return
  154. 3000 rem ----- grafik zeigen ----------
  155. 3010 if w=0 then print"da fehlen noch die werte!":for i=0 to 500:next i:return
  156. 3020 sys 49152:sys 49202,6,0:rem grafik ein
  157. 3021 rem graphic1:return
  158. 3030 get a$:if val(a$)<>4 then 3020
  159. 3040 goto 4010:rem textmodus ein
  160. 4000 rem ----- textmodus ein ----------
  161. 4010 sys 49242:rem textmodus ein
  162. 4011 rem if peek(238)=79 then graphic5:else graphic0
  163. 4020 return
  164. 5000 rem ---- polynom-berechnung ----------------------
  165. 5002 gosub 4000:rem textmodus
  166. 5004 if w=0 then print"da fehlen die werte!":for i=0 to 500:next i:return
  167. 5010 printchr$(147)chr$(17)chr$(17)"welchen grad soll das polynom haben ?"
  168. 5020 print:print"maximal erlaubt ist ein polynom ":print,g".grades ."
  169. 5030 print:print"bei overflow-error ist der":print"wiedereinstieg ins programm"
  170. 5040 print"moeglich mit 'goto 400' !"
  171. 5050 print:input"polynomgrad=";g2:ifg2>gthen5020
  172. 5060 fori=1tog2+2:t(i)=0:a(i)=0:a(abs(2*i-3))=0:fork=1tog+1:r(k,i)=0:nextk:nexti
  173. 5070 a(1)=n:zz=0:m=0:s=0:hh=0:q=0:p=0:a0=0
  174. 5080 fori=1ton:forl=2to2*g2+1:a(l)=a(l)+w(0,i)^(l-1):nextl
  175. 5090 fork=1tog2+1:r(k,g2+2)=t(k)+w(1,i)*w(0,i)^(k-1)
  176. 5100 t(k)=t(k)+w(1,i)*w(0,i)^(k-1):nextk:t(g2+2)=t(g2+2)+w(1,i)^2:nexti
  177. 5110 fori=1tog2+1:fork=1tog2+1:r(i,k)=a(i+k-1):nextk:nexti
  178. 5120 fori=1tog2+1:fork=itog2+1:ifr(k,i)<>0then5150
  179. 5130 gosub 4000:rem textmodus
  180. 5140 print"keine eindeutige loesung":return
  181. 5150 forl=1tog2+2:s=r(i,l):r(i,l)=r(k,l):r(k,l)=s:nextl
  182. 5160 m=1/r(i,i):forl=1tog2+2:r(i,l)=m*r(i,l):nextl
  183. 5170 fork=1tog2+1:ifk=ithen5190
  184. 5180 m=-r(k,i):forl=1tog2+2:r(k,l)=r(k,l)+m*r(i,l):nextl
  185. 5190 nextk:nexti:a0=1:printchr$(147)
  186. 5200 p=0:fori=2tog2+1:p=p+r(i,g2+2)*(t(i)-a(i)*t(1)/n):nexti
  187. 5210 q=t(g2+2)-t(1)^2/n:zz=q-p:b=n-g2-1:hh=p/q:ifb=0thenb=1e-23
  188. 5215 a0=1:gosub 4000:printchr$(147)
  189. 5220 print"das polynom "g2".grades ist:":print:printtab(5)"y=a0+a1*x+a2*x^2+..."
  190. 5230 print:printtab(3)"konstante a0="r(1,g2+2):fori=1tog2
  191. 5240 printtab(3)"koeffizient a"i"="r(i+1,g2+2):nexti:print
  192. 5250 printtab(3)"korrelationskoeffizient=":print,hh
  193. 5260 print:printtab(3)"standardabweichung=":print,sqr(abs(zz/b))
  194. 5270 print:printtab(3)"grafik...taste druecken ! (_ = menue)"
  195. 5280 geta$:if a$=""then 5280
  196. 5290 if a$="_" then return
  197. 5300 r=1:sys 49152:sys 49202,6,0:rem grafik ein
  198. 5301 rem gosub 3000
  199. 5310 for i=lx to mx step (mx-lx)/100
  200. 5320 p=r(1,g2+2):gosub 7010:rem funktionswert berechnen
  201. 5330 x1=fnx(i):y1=fny(p):if y1<0 then 5350
  202. 5340 sys 49266 x1,y1,1:rem punkt zeichnen
  203. 5341 rem draw 1,x1,y1
  204. 5350 next i
  205. 5360 get a$:if a$ ="" then 5360
  206. 5370 if a$="_" then gosub 4000:return
  207. 5380 if r=1 then r=0:gosub 4000:goto5360
  208. 5390 if r=0 then r=1:sys 49152:sys 49202,6,0
  209. 5391 rem if r=0 then r=1:gosub 3000
  210. 5400 goto 5360
  211. 6000 rem ------ werte berechnen ----------
  212. 6010 gosub 4000:rem textmodus
  213. 6020 if a0=0 and w=0 then print"bitte geben sie zuerst werte ein und"
  214. 6030 if a0=0 then print"bitte die kurve berechnen!":for i=0 to 500:next i:return
  215. 6040 print:print" auf der basis der regressionskurve"
  216. 6050 print"koennen beliebige werte berechnet werden"
  217. 6060 print:print" zurueck zum menue kommen sie durch _":print
  218. 6070 input"wert x = ";a$
  219. 6080 if a$ ="_" then return
  220. 6090 i=val(a$)
  221. 6100 p=r(1,g2+2):gosub7010
  222. 6110 printchr$(145),,"y = "p
  223. 6120 goto 6060
  224. 7000 rem ----- polynomwert berechnen -----
  225. 7010 forj=1tog2:p=p+r(j+1,g2+2)*i^j:nextj:return
  226.