10 REM Guido Petz (1098), 1986
20 REM Ickelinj{ra ekvationssystem, ABC80, BASICII-kortet,80 tecken, ABC 802
30 REM ABC80 ger pga.avrundningsfel d}liga resultat vid vissa problem,speciellt
40 REM problem med m}nga variabler!
50 REM Definiera systemet p} rad 2000 som F(i)=fi(x(1),..x(n)) i=1..n
60 REM Programmet ger m|jlighet att anv{nda en grov inb{ddningsteknik dvs
70 REM att l|sa en f|ljd av ekvationssystem som som konvergerar till det
80 REM definierade systemet. Kan vara mycket tidskr{vande !!
90 REM Med cykler menas antalet ggr den LU-faktoriserade iterationsmatrisen
100 REM skall anv{ndas innan den uppdateras och LU-faktoriseras igen.
110 REM Exemplet med 10 variabler ger l|sningar f|r (1,1,1,1,1,1,1,1,1,1 )
120 REM (0,1,0,1,1,0,1,1,1,0) och (-1,1,1,1,-1,0,1,0,1,2) som startvektor!
130 REM Programmet anv{nder en generalisering av Steffensen's metod till N
140 REM dimensioner och de resulterande iterationsmatriserna LU-faktoriseras
150 REM med Crout's metod. Se Bj|rk-Dahlqvist "Numerical methods"!!
160 OPEN 'v24:vsa29c72.1' ASFILE 1
170 ONERRORGOTO 170 : ; CHR$(12)'Ickelinj{ra ekvationssystem!' : ; STRING$(28,45)
180 ; 'Hur m}nga ekvationer:'; : INPUT N% : DIM A(N%,N%),F(N%),B(N%),X(N%),X1(N%),Y(N%),P%(N%)
190 REM parameter f|r ber{kning av iterations matrisen beror p} singel
200 REM eller double precision p} ABC800
210 IF PEEK(-226)=8% E1=1E-14 : G$='1E-14' ELSE E1=.00001 : G$='1E-5'
220 REM Input startvektor+ekvationsevaluering!
230 C%=0% : P1%=0% : F%=0% : ONERRORGOTO 260 : GOSUB 770 : ; 'Printer(j)';
240 GET A$ : ; : IF A$='j' P1%=1% : F%=-1% : GOSUB 740
250 GOTO 230
260 ; 'Inb{ddning(j)'; : GET A$ : ; : IF A$<>'j' 280 ELSE B%=-1% : ; 'Hur m}nga system:'; : INPUT N1%
270 M=1/N1% : FOR I%=1% TO N% : Y(I%)=X(I%) : NEXT I%
280 ONERRORGOTO 290 : I2%=1% : ; 'Noggranhet:('+G$+')'; : INPUT E : GOTO 300
290 E=E1 : ; E : REM default
300 ONERRORGOTO 310 : ; 'Antal iterationer (20)'; : INPUT L% : GOTO 320
310 L%=20% : ; L% : REM default
320 ONERRORGOTO 330 : ; 'Antal cykler (optimerad)'; : INPUT I5% : GOTO 340
330 ONERRORGOTO 0 : C%=-1% : ; 'Optimerad!'
340 ; 'Mellanresultat(j)'; : GET F$ : ; : IF NOT B% 360
350 K=M : FOR I4%=1% TO N1%
360 B=1 : A%=0% : I0%=0% : GOSUB 2000 : IF B% GOSUB 730
370 FOR D%=1% TO L% : REM Ber{kna euklidiska funktionsvektornormen (euklidisk)
380 FOR I%=1% TO N% : B(I%)=F(I%) : X1(I%)=X(I%) : NEXT I% : GOSUB 790
390 G1=G : IF I0%<I2% AND A% 530
400 REM Ber{kna iterationsmatrisen!!
410 I0%=0% : A%=0% : FOR J%=1% TO N% : H=(SGN(B(J%))+.01)*(ABS(B(J%))*.0001+E1)
420 X(J%)=X(J%)+H : GOSUB 2000 : IF B% GOSUB 730
430 FOR I%=1% TO N% : A(I%,J%)=(F(I%)-B(I%))/H : NEXT I% : X(J%)=X(J%)-H : NEXT J%
440 REM Partiell pivotering + LU-faktorisering enl. Crout!
450 FOR K%=1% TO N% : T=0 : FOR I%=K% TO N% : H=A(I%,K%) : FOR P%=1% TO K%-1%
460 H=H-A(I%,P%)*A(P%,K%) : NEXT P% : A(I%,K%)=H : IF ABS(H)>T THEN T=ABS(H) : I1%=I%
470 NEXT I% : P%(K%)=I1% : IF I1%=K% 490
480 FOR J%=1% TO N% : T=A(K%,J%) : A(K%,J%)=A(I1%,J%) : A(I1%,J%)=T : NEXT J%
490 IF A(K%,K%)<>0% T=1/A(K%,K%) ELSE T=E
500 FOR I%=K%+1% TO N% : A(I%,K%)=T*A(I%,K%) : NEXT I% : FOR J%=K%+1% TO N% : H=A(K%,J%)
510 FOR P%=1% TO K%-1% : H=H-A(K%,P%)*A(P%,J%) : NEXT P% : A(K%,J%)=H : NEXT J% : NEXT K% : A%=-1%
520 REM Multiplikation av h|gerled med LU-faktoriserad iterationsmatris!
530 FOR K%=1% TO N% : T=B(P%(K%)) : B(P%(K%))=B(K%) : B(K%)=T : H=B(K%) : FOR P%=1% TO K%-1%
540 H=H-A(K%,P%)*B(P%) : NEXT P% : B(K%)=H : NEXT K% : B=0 : FOR I%=N% TO 1% STEP -1%
550 FOR J%=N% TO I%+1% STEP -1% : B(I%)=B(I%)-A(I%,J%)*B(J%) : NEXT J% : IF A(I%,I%)=0 THEN A(I%,I%)=E
560 ONERRORGOTO 570 : B(I%)=B(I%)/A(I%,I%) : B=B+B(I%)*B(I%)
570 ONERRORGOTO 0 : NEXT I% : B=SQR(B) : A=10 : IF B=0 THEN 670
580 REM Inkrementvektorn divideras med 10 och subtraheras fr}n f|rra
590 REM approximationen tills den nya funktionsvektorns norm {r mindre {n
600 REM den gamlas!
610 A=A*.1 : FOR I%=1% TO N% : X(I%)=X1(I%)-B(I%)*A : NEXT I% : GOSUB 2000 : IF B% GOSUB 730
620 GOSUB 790 : IF G>G1 AND A>1E-15 THEN 610
630 I0%=I0%+1% : IF F$='j' GOSUB 740 : GOTO 660 ELSE ; CHR$(12)CUR(10,6); : IF B% ; 'System:'I4%;
640 ; ' Iteration:'D%' av'L%' Fel:'B
650 REM Antal cykler {r 7 om inkrementvektorns norm B < 1 annars 1
660 IF I0%<I2% 670 ELSE IF C% I2%=-(7*(B<1)+(B>=1)) ELSE I2%=I5%
670 IF B<E THEN 680 ELSE NEXT D%
680 IF B% K=K+M : NEXT I4% : REM inbeddning n{sta system!
690 I4%=N1% : IF F$<>'j' OR P1%=1% ; #P1%CHR$(12) : GOSUB 740 : ; CHR$(7)
700 IF P1%=0% ; 'Printer(j)'; : GET A$ : ; : IF A$='j' P1%=1% : GOTO 690
710 P1%=0% : ; 'Forts{tta(j)'; : GET A$ : ; : B%=0% : IF A$='j' 230
720 END
730 FOR I%=1% TO N% : F(I%)=F(I%)*K+(1-K)*(X(I%)-Y(I%)) : NEXT I% : RETURN
740 IF B% ; #P1%'System:'I4%;
750 IF F% ; #P1%' Start:' ELSE ; #P1%'Iteration:'D%
760 FOR I%=1% TO N% : ; #P1%'X'I%'='X(I%);TAB(35%)'F'I%'='F(I%) : NEXT I% : RETURN
770 FOR I%=1% TO N% : ; 'X'I%'='; : INPUT C$ : X(I%)=VAL(C$) : NEXT I%
780 GOSUB 2000 : ; : FOR I%=1% TO N% : ; 'F'I%'='F(I%) : NEXT I% : RETURN
790 ONERRORGOTO 800 : G=0 : FOR I%=1% TO N% : G=G+F(I%)*F(I%) : NEXT I% : G=SQR(G)
800 RETURN
2000 ONERRORGOTO 2110
2010 F(1)=X(1)*ATN(X(2)-X(10))+COS(X(7)*X(8))-X(3)*X(4)-TAN(X(8))+2*X(1)
2020 F(2)=X(2)*X(2)*EXP(X(5)*X(6))+EXP(-X(8)*X(9))-X(1)*X(7)-X(3)-X(9)
2030 F(3)=SIN(1-X(2)*X(10))+X(2)*X(7)*X(3)-X(1)*X(1)*X(10)*X(10)+TAN(X(3))-TAN(X(8))
2040 F(4)=X(2)*X(10)-X(8)*X(9)+SIN(X(5))-X(4)*X(7)*X(7)
2050 F(5)=X(5)*EXP(X(8))-X(4)*SIN(X(7))-COS(X(10))-X(5)
2060 F(6)=EXP(COS(X(5)))-EXP(X(3))-X(1)*X(2)*X(10)+X(6)*X(6)+X(1)
2070 F(7)=X(6)*X(7)+X(2)*X(9)-X(1)*SIN(X(8))-X(4)*SIN(X(7))
2080 F(8)=ATN(1-X(9))-COS(X(3)-X(6))-X(2)*X(5)+2*X(8)
2090 F(9)=X(9)*EXP(X(4)*X(6))-TAN(X(2)*X(5))-X(2)*X(10)+X(1)*X(9)-X(6)
2100 F(10)=X(3)*X(8)*X(10)-X(4)*X(7)*X(10)-X(1)*X(10)+X(2)*X(9)
2110 ONERRORGOTO 0 : RETURN