Gleichungen vom Typ
F[1]=f1(x[1],x[2],x[3]...,
10 PROGRAM FEHLBERG;
20
30 {$L-}
40
50 CONST GRAD=3;
60 N1=4;
70
80 TYPE
90
100 VEKTOR=ARRAY[1..N1] OF REAL;
110
120 VAR
130
140 F,XO,P :VEKTOR;
150
160 K1,K2,K3,K4,K5,K6:VEKTOR;
170
180 EPS,TA,TE,TS :REAL;
190
200 N: INTEGER;
210
220 A :ARRAY[1..N1,1..500] OF REAL;
230
240 CR,B :CHAR;
250
260
270 PROCEDURE TEXT;
280
290 VAR I:INTEGER;
300
310 BEGIN
320 PAGE;WRITELN;WRITELN;
330 WRITELN(' Loesung eines DGLS');
340 WRITELN;
350 WRITELN(' Gleichungen vom Typ');
360 WRITELN;
370 WRITELN(' F[1]=f1(x[1],x[2],x[3]...,T)');
380 WRITELN;
390 WRITELN(' F[2]=f2(x[1],x[2],x[3]...,T)');
400 WRITELN;
410 WRITELN(' F[3]=f3(x[1],x[2],x[3]...,T)');
420 WRITELN;WRITELN;
430 WRITELN(' Eingabe des DGLS ab Zeile 1140');
440 WRITELN;
450 WRITE(' --> ENTER');
460 READ(CR);
470
480 END; {PROC. TEXT}
490
500
510 PROCEDURE ZEIT;
520
530 BEGIN
540 PAGE;WRITELN;WRITELN;
550 WRITE(' Intervallanfang TA=');READ(TA);
560
570 WRITELN;
580 WRITE(' Intervallende TE=');READ(TE);
590 WRITELN;
600 WRITE(' Druckschrittweite TS=');READ(TS);
610 WRITELN;
620 WRITE(' max.lokaler Fehler e=');READ(EPS);
630
640 END;{PROC.ZEIT}
650
660
670
680
690 PROCEDURE ANFANG;
700
710 VAR I:INTEGER;
720
730 BEGIN
740 PAGE;WRITELN;WRITELN;
750 WRITELN('Eingabe des Startvektors !');
760 WRITELN;WRITELN;
770 FOR I:=1 TO GRAD DO
780 BEGIN
790 WRITE(' Xo[',I,']='); READ(XO[I]);
800 END;
810 END; {PROC.ANFANG}
820
830
840
850 PROCEDURE PARAMETER;
860
870
880 CONST PA=2;
890
900
910 VAR I :INTEGER;
920
930 BEGIN
940 PAGE;WRITELN;WRITELN;
950 WRITELN('alte Parameter:');
960 FOR I:=1 TO PA DO
970 WRITELN(' P',I,'=',P[I]:7:4);
980 WRITELN;WRITELN;
990 WRITELN('neue Parameter:');
1000 FOR I:=1 TO PA DO
1010 BEGIN
1020 WRITE(' P',I,'= ');
1030 READ(P[I]);
1040 END;
1050
1060 END; {PROC.PARAMETER}
1070
1080
1090 PROCEDURE DGLS(X:VEKTOR;T:REAL;VAR F:VEKTOR);
1100
1110
1120 CONST AL=0.687;
1130 BE=9.85;
1140 L1=0.025;
1150 L2=0.075;
1160
1170
1180
1190 VAR Q9,Q8 :REAL;
1200
1210
1220 BEGIN
1230
1240 IF ABS(X[1])<1E-4 THEN
1250 BEGIN
1260 X[1]:=0;
1270 Q8:=0;
1280 Q9:=0;
1290 END
1300 ELSE
1310 BEGIN
1320 Q8:=SQR(SQR(SQR(X[1])));
1330 Q9:=Q8*X[1];
1340 END;
1350
1360 F[1]:=X[2];
1370
1380 F[2]:=X[3];
1390
1400 F[3]:=-P[1]*X[3]-AL*(1+ BE*Q8)*X[2]-L1*X[1]-L2* Q9+P[2]*SIN(T);
1410
1420 END; {PROC. DGLS}
1430
1440
1450
1460
1470
1480 PROCEDURE STEP(T0,H:REAL;X0:VEKTOR;VAR X1:VEKTOR; VAR MAXF:REAL);
1490
1500 VAR
1510 I:INTEGER;
1520 T:REAL;
1530 DEL,X:VEKTOR;
1540
1550 BEGIN
1560
1570 X:=X0;
1580 T:=T0;
1590 DGLS(X,T,K1);
1600 FOR I:=1 TO GRAD DO
1610 BEGIN
1620 K1[I]:=K1[I]*H;
1630 X[I]:=X0[I]+2*K1[I]/9;
1640 END;
1650 T:=T0+2*H/9;
1660 DGLS(X,T,K2);
1670 FOR I:=1 TO GRAD DO
1680 BEGIN
1690 K2[I]:=K2[I]*H;
1700 X[I]:=X0[I]+K1[I]/12+K2[I]/4;
1710 END;
1720 T:=T0+H/3;
1730 DGLS(X,T,K3);
1740 FOR I:=1 TO GRAD DO
1750 BEGIN
1760 K3[I]:=K3[I]*H;
1770 X[I]:=X0[I]+(69*K1[I]-243*K2[I]+270*K3[I])/128;
1780 END;
1790 T:=T0+0.75*H;
1800 DGLS(X,T,K4);
1810 FOR I:=1 TO GRAD DO
1820 BEGIN
1830 K4[I]:=K4[I]*H;
1840 X[I]:=X0[I]+(-85*K1[I]+405*K2[I]-324*K3[I]+64*K4[I])/60;
1850 END;
1860 T:=T0+H;
1870 DGLS(X,T,K5);
1880 FOR I:=1 TO GRAD DO
1890 BEGIN
1900 K5[I]:=K5[I]*H;
1910 X[I]:=X0[I]+(65*K1[I]-135*K2[I]+351*K3[I]+64*K4[I]+15*K5[I])/432;
1920 END;
1930 T:=T0+5*H/6;
1940 DGLS(X,T,K6);
1950 MAXF:=0;
1960 FOR I:=1 TO GRAD DO
1970 BEGIN
1980 K6[I]:=K6[I]*H;
1990 X1[I]:=X0[I]+(20*K1[I]+81*K3[I]+64*K4[I]+15*K5[I])/180;
2000 DEL[I]:=ABS((2*K1[I]-3*K3[I]+64*K4[I]+15*K5[I]-72*K6[I])/300);
2010 IF MAXF<DEL[I] THEN MAXF:=DEL[I];
2020 END;
2030
2040 END; {PROC.STEP}
2050
2060
2070
2080
2090
2100 PROCEDURE WAHL(VAR XA,XB:INTEGER);
2110
2120 BEGIN
2130
2140 PAGE;
2150 WRITELN;WRITELN;
2160 WRITELN('Welche Groessen zur Ausgabe ?');
2170 WRITELN;
2180 WRITELN;
2190 WRITELN('X[1] --> 1');
2200 WRITELN;
2210 WRITELN('X[2] --> 2');
2220 WRITELN;
2230 WRITELN('X[n] --> n');
2240 WRITELN;
2250 WRITELN(' T --> ',N1);
2260 WRITELN;WRITELN;
2270 WRITE(' X-Achse :=');READ(XA);
2280 WRITELN;WRITELN;
2290 WRITE(' Y-Achse :=');READ(XB);
2300
2310 END; {PROC. WAHL}
2320
2330 PROCEDURE AUSGABE;
2340
2350 VAR I,J,K,ZE :INTEGER;
2360
2370 BEGIN
2380
2390 PAGE;
2400 WRITELN;
2410 WRITELN(' numerische Ausgabe');
2420 WRITELN(' ******************');
2430 WRITELN;
2440 WRITELN('Welche Zustandsvariablen:');
2450 WRITELN;
2460 WRITE('Xi --> i:=');
2470 READ(K);
2480 WRITELN;WRITELN;
2490 WRITE('Xj --> j:=');
2500 READ(J);PAGE;
2510 ZE:=0;
2520 WRITELN(' X',K,' X',J,' Zeit');
2530 WRITELN;
2540 FOR I:=1 TO N DO
2550 BEGIN
2560 ZE:=ZE+1;
2570 WRITELN(A[K,I]:10:5,A[J,I]:10:5,A[N1,I]:10:2);
2580 WRITELN;
2590 IF ZE=13 THEN
2600 BEGIN
2610 WRITE('--> ENTER'); READ(CR);
2620 ZE:=0;
2630 PAGE;
2640 WRITELN(' X',K,' X',J,' Zeit');
2650 WRITELN;
2660 END;
2670 END;
2680 WRITE('--> ENTER');
2690 READ(CR);
2700
2701
2702
2703
2704
2710 END; {PROC. AUSGABE}
2720
2730 PROCEDURE GRAPH;
2740
2750 VAR
2760 T1,I,J,K,D,M :INTEGER;
2770 Y1,Y2,Z1,Z2 :REAL;
2780 F1,P,R,S :REAL;
2790 XA,XB :INTEGER;
2800
2810
2820 BEGIN
2830 WAHL(XA,XB);
2840 Y1:=1E30;Y2:=-1E30;
2850 Z1:=1E30;Z2:=-1E30;
2860 FOR I:=1 TO N DO
2870 BEGIN
2880 IF A[XA,I]<Y1 THEN Y1:= A[XA,I];
2890 IF A[XA,I]>Y2 THEN Y2:= A[XA,I];
2900 IF A[XB,I]<Z1 THEN Z1:= A[XB,I];
2910 IF A[XB,I]>Z2 THEN Z2:= A[XB,I];
2920 END;
2930 OUT(4,'8');
2940 PAGE;
2950 P:=ABS(Y2-Y1);
2960 R:=ABS(Z2-Z1);
2970 P:=124/P;
2980 R:=124/R;
2990 IF XA<N1 THEN
3000 BEGIN
3010 IF P<R THEN R:=P;
3020 IF R<P THEN P:=R;
3030 END;
3040 FOR K:=1 TO N DO
3050 BEGIN
3060 I:=ENTIER((A[XA,K]-Y1)*P);
3070 M:=ENTIER(I/4);
3080 F1:=(A[XB,K]-Z1)*R;
3090 D:=ENTIER(F1/4);
3100 J:=220+I-4*M-4*ENTIER(F1)+16*D;
3110 POKE(-4128-32*D+M,CHR(J));
3120 END;
3130 F1:=-Z1*R;
3140 D:=ENTIER(F1/4);
3150 I:=ENTIER(-Y1*P);
3160 M:=ENTIER(I/4);
3170 FOR T1:=1 TO 32 DO
3180 BEGIN
3190 POKE(-4128-32*T1+M,CHR(161));
3200 POKE(-4129-32*D+T1,CHR(160));
3210 END;
3220 WRITE('Ent.');READ(CR);
3230 END; {PROC. GRAPH}
3240
3250
3260
3270
3280
3290 PROCEDURE INTEGRAL;
3300
3310 VAR
3320 H,H1,MAXF:REAL;
3325 T,TV :REAL;
3330 J2 :INTEGER;
3340 X0,X1 :VEKTOR;
3350
3360 BEGIN
3370 XO[N1]:=TA;
3380 T:=TA;
3390 FOR J2:=1 TO N1 DO
3400 A[J2,1]:=XO[J2];
3410 X0:=XO;
3420 N:=1;
3430 TV:=T+TS;
3440 H:=TS*0.5;
3460 REPEAT
3470 REPEAT
3480 STEP(T,H,X0,X1,MAXF);
3481
3490 IF MAXF>EPS THEN H:=H *0.5;
3500 UNTIL MAXF<=EPS;
3510 IF T+H>TV THEN
3520 BEGIN
3530 H1:=TV-T;
3540 STEP(T,H1,X0,X1,MAXF);
3550 X0:=X1;N:=N+1;
3560 T:=TV;TV:=T+TS;
3570 PAGE;WRITE('T=',T:5:2,' N=',N);
3580 FOR J2:=1 TO GRAD DO
3590 A[J2,N]:=X1[J2];
3600 A[N1,N]:=T;
3610 END
3620 ELSE
3630 BEGIN
3640 T:=T+H;
3650 X0:=X1;
3660 END;
3670 IF MAXF<(EPS/32) THEN H:=H*1.5;
3680 UNTIL T+TS>=TE;
3690 OUT(4,'8');
3700
3710 END; {PROC.INTEGRAL}
3720
3730
3740 {HAUPTPROGRAMM}
3750
3760 BEGIN
3770 TEXT;
3780 REPEAT
3790 PAGE;WRITELN;
3800 WRITELN(' Menue');
3810 WRITELN(' -----');
3820 WRITELN;WRITELN;
3830 WRITELN('Zeitintervall -->Z');
3840 WRITELN;
3850 WRITELN('Parameterwahl -->P');
3860 WRITELN;
3870 WRITELN('Startvektor -->S');
3880 WRITELN;
3890 WRITELN('Integration -->I');
3900 WRITELN;
3910 WRITELN('graph.Ausgabe -->G');
3920 WRITELN;
3930 WRITELN('numer.Ausgabe -->N');
3940 WRITELN;
3950 WRITELN('Ende -->E');
3960 WRITELN;
3970 WRITE(' skip Key -->');
3980 B:=INCH;
3990 CASE B OF
4000
4010 'Z': ZEIT;
4020 'S': ANFANG;
4030 'P': PARAMETER;
4040 'I': INTEGRAL;
4050 'G': GRAPH;
4060 'N': AUSGABE
4070 END;
4080
4090 UNTIL B='E';
4100
4110 END.