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 >
Wrap
Commodore BASIC
|
2022-10-26
|
8KB
|
226 lines
1 rem ----- nachladen grafik -------
2 if a$="n" then 100
3 printchr$(147)chr$(17)"grafik nachladen (j/n)"
4 get a$:if a$="" then 4
5 if a$="n" then 100
6 a$="n":load"sysgraf.obj",8,1
10 rem*********************************
20 rem* *
30 rem* kurvenanpassung *
40 rem* *
50 rem* statistik-programm zur *
60 rem* regressionsanalyse mittels *
65 rem* polynomen incl.grafik *
70 rem* heimo ponnath hamburg 1987 *
80 rem* c64 - version *
90 rem*********************************
100 clr:sys 49152:sys 49242:rem grafikspeicher sichern
105 rem graphic1,1:gosub 4000
110 rem ----- variable ----------------
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
130 bb=319:bh=199:w=0:mx=-1e12:lx=1e12:my=-1e12:ly=1e12:r=0
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
150 a$="":b$=""
160 rem
170 deffn g(n)=int(.2*sqr(5*(a-10*n-17))-4):rem max. polynomgrad berechnen
180 deffn x(x)=ra*x+ta:deffn y(y)=rd*y+tb
190 rem
200 rem ----- titel,erklaerung --------
210 poke 53280,0:poke 53281,0:printchr$(30)
220 printchr$(147)chr$(18)" polynomanpassung "chr$(146)
230 print
240 print" durch eine anzahl von n punkten aus"
250 print"wertepaaren legt dieses programm das am"
260 print"besten angepasste polynom der form"
270 print"y=a0+a1*x+a2*x^2+a3*x^3+... . der grad"
275 print"des polynoms ist frei waehlbar. der"
277 print"korrelationskoeffizient r und die"
280 print"standardabweichung s werden angegeben"
290 print"und sie koennen beliebige y-werte"
300 print"aus eingegebenen x-werten berechnen.":print
310 print" ein scatterdiagramm und die ermittelte"
320 print"kurve werden gezeichnet. auf diese"
330 print"weise kann die qualitaet der anpassung"
340 print"eingeschaetzt werden.":print
350 print" sogenannte ausreisser-werte sollten"
360 print"vor einer genaueren berechnung noch"
370 print"entfernt werden.":print:print
380 printchr$(18)"taste druecken!"chr$(146)
390 get a$:if a$="" then 390
400 rem ----- hauptmenue --------------
410 printchr$(147):print:print:print:print
420 printtab(4)"werte von hand eingeben.......1":print
430 printtab(4)"werte aus datei lesen.........2":print
440 printtab(4)"grafik zeigen.................3":print
450 printtab(4)"textmodus einschalten.........4":print
460 printtab(4)"polynomfunktion berechnen.....5":print
470 printtab(4)"werte berechnen...............6":print
480 printtab(4)"programmende..................7":print:print
490 printtab(10)chr$(18)"bitte waehlen sie!"chr$(146)
500 get a$:if val(a$)<1 or val(a$)>7 then 500
510 printchr$(147):if val(a$)=7 then end
520 on val(a$) gosub 1000,2000,3000,4000,5000,6000
530 goto 410
540 rem ----- ende hauptprogramm ------
1000 rem ----- werte von hand ---------
1005 gosub 4000:rem textmodus
1010 if w=1 then print"werte schon vorhanden!":for i=0 to 500:next i:return
1020 w=1
1030 print"wieviele werte werden verwendet ?":inputn:print
1040 dim w(1,n)
1043 a=fre(0)-2000:rem freier speicherplatz c64
1045 rem a=fre(1)-2000
1047 g=fng(n):rem maximaler polynomgrad
1050 print"bitte wertepaare eingeben!":print
1060 for i=1 to n
1070 printi,"x=";:inputw(0,i):printchr$(145),,"y=";:inputw(1,i):print
1080 gosub 1300:rem zwischenwerte berechnen
1090 next i
1100 printchr$(147):print:print"sollen die werte gespeichert werden?"
1110 get a$:if a$<>"j" and a$<>"n" then 1110
1120 if a$="n" then 1190
1130 print:print"name der datei (11 zeichen)";:input b$
1140 b$=left$(b$,11)+".dat"+",s,w"
1150 open 1,8,2,b$
1160 print#1,n
1170 for i=1 to n:print#1,w(0,i):print#1,w(1,i):next i
1180 close 1
1190 gosub 1500:rem scatterdiagramm zeichnen
1200 return
1300 rem --- zwischenwerte berechnen --
1360 if w(0,i)>mx then mx=w(0,i)
1370 if w(0,i)<lx then lx=w(0,i)
1380 if w(1,i)>my then my=w(1,i)
1390 if w(1,i)<ly then ly=w(1,i)
1400 return
1500 rem --- scatterdiagramm ----------
1505 dim a(2*g+1),r(g+1,g+2),t(g+2):rem arrays fuer berechnungen
1510 for i=1 to n-1:rem sortieren nach x
1520 for j=i+1 to n
1530 if w(0,i)<=w(0,j) then 1560
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)
1550 w(0,j)=w(0,0):w(1,j)=w(1,0)
1560 next j:next i
1570 sys 49152:sys 49180:sys 49202,6,0:rem grafik loeschen farbe
1580 rem graphic1,1:color0,1:color1,7
1590 sys 49352,0,0,319,0,1:sys 49352,319,0,319,199,1
1600 rem draw 1,0,0 to 319,0 to 319,199 to 0,199 to 0,0
1610 sys 49352,319,199,0,199,1:sys 49352,0,199,0,0,1:rem rahmen
1620 xu=lx-(mx-lx)*.02:xo=mx+(mx-lx)*.02
1630 yu=ly-(my-ly)*.02:yo=my+(my-ly)*.02
1640 ra=bb/(xo-xu):rd=-bh/(yo-yu)
1650 ta=-bb*xu/(xo-xu):tb=bh*yo/(yo-yu)
1660 for i=1 to n
1670 x=fnx(w(0,i)):y=fny(w(1,i))
1680 sys49352,x-3,y,x+3,y,1:sys49352,x,y-3,x,y+3,1:rem kreuz
1681 rem draw1,x-3,y to x+3,y:draw1,x,y-3 to x,y+3
1690 next i
1700 get a$:if a$="" then 1700
1710 sys 49242:rem textmodus
1711 rem if peek(238)=79 then graphic5:else graphic0
1720 print"xu = "lx,"xo = "mx"
1730 [153]"yu = "ly,"yo = "my"
1740 get a$:if a$="" then 1740
1750 return
2000 rem ----- werte aus datei --------
2005 gosub 4000:rem textmodus
2010 if w=1 then print"werte schon vorhanden!":for i=0 to 500:next i:return
2020 w=1
2030 print" die datei muss ein bestimmtes format"
2040 print"haben: 1.anzahl der wertepaare"
2050 print" 1.wert x, 1.wert y"
2060 print" 2.wert x, 2.wert y ...":print
2070 print"diese dateien werden unter menuepunkt 1"
2080 print"erstellt. sie tragen die endung .dat .":print
2090 print" alles klar..1 ach soo..2"
2100 get a$:if val(a$)<1 or val(a$)>2 then 2100
2110 if val(a$)=2 then w=0:return
2120 print:print"wie heisst denn die datei (endung .dat)"
2130 input b$
2140 b$=b$+",s,r"
2150 open1,8,2,b$
2160 input#1,n
2170 dim w(1,n)
2173 a=fre(0)-2000:rem freier speicherplatz c64
2175 rem a=fre(1)-2000
2177 g=fng(n):rem maximaler polynomgrad
2180 for i=1 to n
2190 input#1,w(0,i):input#1,w(1,i)
2200 gosub 1300:rem zwischenwerte berechnen
2210 next i
2220 close 1
2230 gosub 1500:rem scatterdiagramm
2240 return
3000 rem ----- grafik zeigen ----------
3010 if w=0 then print"da fehlen noch die werte!":for i=0 to 500:next i:return
3020 sys 49152:sys 49202,6,0:rem grafik ein
3021 rem graphic1:return
3030 get a$:if val(a$)<>4 then 3020
3040 goto 4010:rem textmodus ein
4000 rem ----- textmodus ein ----------
4010 sys 49242:rem textmodus ein
4011 rem if peek(238)=79 then graphic5:else graphic0
4020 return
5000 rem ---- polynom-berechnung ----------------------
5002 gosub 4000:rem textmodus
5004 if w=0 then print"da fehlen die werte!":for i=0 to 500:next i:return
5010 printchr$(147)chr$(17)chr$(17)"welchen grad soll das polynom haben ?"
5020 print:print"maximal erlaubt ist ein polynom ":print,g".grades ."
5030 print:print"bei overflow-error ist der":print"wiedereinstieg ins programm"
5040 print"moeglich mit 'goto 400' !"
5050 print:input"polynomgrad=";g2:ifg2>gthen5020
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
5070 a(1)=n:zz=0:m=0:s=0:hh=0:q=0:p=0:a0=0
5080 fori=1ton:forl=2to2*g2+1:a(l)=a(l)+w(0,i)^(l-1):nextl
5090 fork=1tog2+1:r(k,g2+2)=t(k)+w(1,i)*w(0,i)^(k-1)
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
5110 fori=1tog2+1:fork=1tog2+1:r(i,k)=a(i+k-1):nextk:nexti
5120 fori=1tog2+1:fork=itog2+1:ifr(k,i)<>0then5150
5130 gosub 4000:rem textmodus
5140 print"keine eindeutige loesung":return
5150 forl=1tog2+2:s=r(i,l):r(i,l)=r(k,l):r(k,l)=s:nextl
5160 m=1/r(i,i):forl=1tog2+2:r(i,l)=m*r(i,l):nextl
5170 fork=1tog2+1:ifk=ithen5190
5180 m=-r(k,i):forl=1tog2+2:r(k,l)=r(k,l)+m*r(i,l):nextl
5190 nextk:nexti:a0=1:printchr$(147)
5200 p=0:fori=2tog2+1:p=p+r(i,g2+2)*(t(i)-a(i)*t(1)/n):nexti
5210 q=t(g2+2)-t(1)^2/n:zz=q-p:b=n-g2-1:hh=p/q:ifb=0thenb=1e-23
5215 a0=1:gosub 4000:printchr$(147)
5220 print"das polynom "g2".grades ist:":print:printtab(5)"y=a0+a1*x+a2*x^2+..."
5230 print:printtab(3)"konstante a0="r(1,g2+2):fori=1tog2
5240 printtab(3)"koeffizient a"i"="r(i+1,g2+2):nexti:print
5250 printtab(3)"korrelationskoeffizient=":print,hh
5260 print:printtab(3)"standardabweichung=":print,sqr(abs(zz/b))
5270 print:printtab(3)"grafik...taste druecken ! (_ = menue)"
5280 geta$:if a$=""then 5280
5290 if a$="_" then return
5300 r=1:sys 49152:sys 49202,6,0:rem grafik ein
5301 rem gosub 3000
5310 for i=lx to mx step (mx-lx)/100
5320 p=r(1,g2+2):gosub 7010:rem funktionswert berechnen
5330 x1=fnx(i):y1=fny(p):if y1<0 then 5350
5340 sys 49266 x1,y1,1:rem punkt zeichnen
5341 rem draw 1,x1,y1
5350 next i
5360 get a$:if a$ ="" then 5360
5370 if a$="_" then gosub 4000:return
5380 if r=1 then r=0:gosub 4000:goto5360
5390 if r=0 then r=1:sys 49152:sys 49202,6,0
5391 rem if r=0 then r=1:gosub 3000
5400 goto 5360
6000 rem ------ werte berechnen ----------
6010 gosub 4000:rem textmodus
6020 if a0=0 and w=0 then print"bitte geben sie zuerst werte ein und"
6030 if a0=0 then print"bitte die kurve berechnen!":for i=0 to 500:next i:return
6040 print:print" auf der basis der regressionskurve"
6050 print"koennen beliebige werte berechnet werden"
6060 print:print" zurueck zum menue kommen sie durch _":print
6070 input"wert x = ";a$
6080 if a$ ="_" then return
6090 i=val(a$)
6100 p=r(1,g2+2):gosub7010
6110 printchr$(145),,"y = "p
6120 goto 6060
7000 rem ----- polynomwert berechnen -----
7010 forj=1tog2:p=p+r(j+1,g2+2)*i^j:nextj:return