4823 6164 0000 P ... FEHLBERG
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.