home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / 1987 / 05 / dgl.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1987-04-15  |  4.6 KB  |  137 lines
  1. {---------------------------------------------------------------------------}
  2. {                               DGL.PAS                                     }
  3. {      Loesen von Differentialgleichungssystemen nach Runge-Kutta           }
  4. {---------------------------------------------------------------------------}
  5.  
  6. PROGRAM DGL_Solver;
  7.  
  8. CONST  nm = 5;                                            { Max. DGL-Anzahl }
  9.  
  10. TYPE   Vektor = ARRAY[1..nm] OF REAL;
  11.  
  12. VAR    Y, Yo, dY :Vektor;
  13.        Xstart, Xend, x, dx, h, AbsErr, ymin, ymax :REAL;
  14.        n, i, nmax:INTEGER;
  15.  
  16. {---------------------------------------------------------------------------}
  17. {                  das DGL-System ist hier zu definieren:                   }
  18.  
  19. PROCEDURE DGL (VAR dY, Y: Vektor; x: REAL);
  20.  
  21. CONST  k1=0.04;
  22.        k2=0.03;
  23.  
  24. BEGIN                                          {     n DGLn der Form:       }
  25.   dY[1] := -k1*Y[1];                           { dY[i] := f(x,Y[1],..,Y[n]) }
  26.   dY[2] :=  k1*Y[1] - k2*Y[2];
  27.   dY[3] :=            k2*Y[2];                              { Folgereaktion }
  28. END;
  29.  
  30. {---------------------------------------------------------------------------}
  31.  
  32. PROCEDURE Ausgabe (VAR Y: Vektor; x, h: REAL; n, nmax: INTEGER);
  33.  
  34. CONST xpos=60; ypos=10;                                { Bildschirmposition }
  35. VAR   i: INTEGER;
  36.  
  37. BEGIN
  38.   FOR i:=1 TO n DO
  39.   BEGIN
  40.     GotoXY(xpos,ypos+i); WriteLn('Y[',i,'] =',y[i]);
  41.   END;
  42.   GotoXY(xpos,ypos+n+1); WriteLn('x    =',x);
  43.   GotoXY(xpos,ypos+n+2); WriteLn('h    =',h);
  44.   GotoXY(xpos,ypos+n+3); WriteLn(nmax,' Teilintervalle');
  45. END;
  46.  
  47. {---------------------------------------------------------------------------}
  48.  
  49. PROCEDURE Plot (x, y, xmin, xmax, ymin, ymax: REAL; ch:CHAR);
  50.  
  51. CONST umax=78; vmax=28;                           { Bildschirmbreite/-hoehe }
  52. VAR   u, v: INTEGER;
  53.  
  54. BEGIN
  55.   u :=      Round((x-xmin)/(xmax-xmin)*umax);
  56.   v := vmax-Round((y-ymin)/(ymax-ymin)*vmax);
  57.   GotoXY(u,v); Write(ch);
  58. END;
  59.  
  60. {---------------------------------------------------------------------------}
  61. {                      das Runge-Kutta-Verfahren:                           }
  62.  
  63. PROCEDURE RungeKutta (VAR Y, Yo:Vektor; Xstart, Xend: REAL;
  64.                       VAR h: REAL; AbsErr: REAL;
  65.                           n: INTEGER; VAR nmax: INTEGER);
  66.  
  67. VAR C1, C2, C3, C4, z, dY: Vektor;
  68.     x, xo, yalt, Tol: REAL;
  69.     j, k: INTEGER;
  70.  
  71. BEGIN
  72.   yalt := 0; nmax := 4;
  73.   REPEAT
  74.     h := ABS((Xend-Xstart)/nmax);                            { Schrittweite }
  75.     FOR k := 1 TO n DO Y[k] := Yo[k];                  { Anfangsbedingungen }
  76.     FOR j := 0 TO Pred(nmax) DO                     { Runge-Kutta-Iteration }
  77.     BEGIN
  78.       xo := Xstart+j*h;
  79.       DGL(dY, Y, xo);
  80.       FOR k := 1 TO n DO
  81.       BEGIN
  82.         C1[k] := h*dY[k];
  83.         z[k] := Y[k]+C1[k]/2;
  84.       END;
  85.       DGL(dY, z, (xo+h/2));
  86.       FOR k := 1 TO n DO
  87.       BEGIN
  88.         C2[k] := h*dY[k];
  89.         z[k] := Y[k]+C2[k]/2;
  90.       END;
  91.       DGL(dY, z, (xo+h/2));
  92.       FOR k := 1 TO n DO
  93.       BEGIN
  94.         C3[k] := h*dY[k];
  95.         z[k] := Y[k]+C3[k];
  96.       END;
  97.       DGL(dY, z, (xo+h));
  98.       FOR k := 1 TO n DO
  99.       BEGIN
  100.         C4[k] := h*dY[k];
  101.         y[k] := y[k]+(C1[k]+2*C2[k]+2*C3[k]+C4[k])/6;
  102.       END;
  103.     END;
  104.     Tol := ABS(yalt-y[1]);                               { Abbruchkriterium }
  105.     nmax := nmax*2;                                        { Teilintervalle }
  106.     yalt := Y[1];
  107.   UNTIL Tol <= AbsErr;
  108. END;
  109.  
  110. {---------------------------------------------------------------------------}
  111.  
  112. BEGIN { DGL_Solver }
  113.                                                       { Anfangsbedingungen: }
  114.   n     := 3;                                               { Zahl der DGLn }
  115.   Yo[1] := 2.5;                            { n Anfangsbedingungen der Form: }
  116.   Yo[2] := 0.0;                            {        Yo[i] := const          }
  117.   Yo[3] := 0.0;
  118.   Xstart:= 0.1;                                 { untere Integrationsgrenze }
  119.   Xend  := 120.0;                               { obere  Integrationsgrenze }
  120.   dx    := (Xend-Xstart)/80;                               { x-Schrittweite }
  121.   AbsErr:= 1.0E-3;                                       { Absoluter Fehler }
  122.   ymin  := 0.0;                                  { y-Ausdehnung der Graphik }
  123.   ymax  := 2.5;
  124.  
  125.   ClrScr;
  126.   FOR i := 1 TO n DO Y[i] := Yo[i];
  127.   x := Xstart;
  128.   WHILE x <= Xend DO
  129.   BEGIN
  130.     FOR i := 1 TO n DO Plot(x, Y[i], Xstart, Xend, ymin, ymax, Chr(48+i));
  131.     RungeKutta(Y, Yo, Xstart, x, h, AbsErr, n, nmax);
  132.     Ausgabe(Y, x, h, n, nmax);
  133.     x := x + dx;
  134.   END;
  135.   GotoXY(0,0);
  136. END.
  137.