Turbo Toolbox/Turbo_Toolbox.iso/1987/05/partdgl.pas

Pascal/Delphi Source File | 1987-04-15 | 5.4 KB | 156 lines

{---------------------------------------------------------------------------}
{                            PARTDGL.PAS                                    }
{         Loesen von partiellen DGL-Systemen nach Runge-Kutta               }
{---------------------------------------------------------------------------}

PROGRAM Partial_Dgl_Solver;

CONST  nm = 20;                                    { Max. Stuetzstellenzahl }

TYPE   Vektor = ARRAY[1..nm] OF REAL;

VAR    Y, Yo, dY :Vektor;
       Xstart, Xend, x, dx, h, AbsErr, t, t0, tend, dt, ymin, ymax: REAL;
       n, i, k, nmax: INTEGER;

{****************************************************************************
 *                                                                          *
 *            Schreibweise einer Partiellen DGL in Prozedur DGL             *
 *                                                                          *
 *            dy         d"y                                                *
 *            --   =  k  ---   =  k ( y[i-1] - 2*y[i] + y[i+1] )            *
 *            dt         dx"                                                *
 *                                                                          *
 ****************************************************************************}

PROCEDURE DGL (VAR dY, Y: Vektor; x: REAL);

CONST k=0.0023;
VAR   i: INTEGER;

BEGIN                                             { Partielle DGL der Form: }
  FOR i := 2 TO Pred(n) DO                        { dY[i]:=f(x,Y[1],..,Y[n] }
    dY[i] :=  k * (Y[i-1] - 2*Y[i] + Y[i+1]);           { FOURIER-Gleichung }
  dY[1] :=  Yo[1];
  dY[n] :=  Yo[n];
END;

{---------------------------------------------------------------------------}

PROCEDURE Ausgabe (VAR Y: Vektor; x, h: REAL; n, nmax: INTEGER);

CONST xpos=0; ypos=0;                                  { Bildschirmposition }
VAR   i: INTEGER;

BEGIN
  FOR i := 1 TO n DO
  BEGIN
    GotoXY(xpos,ypos+i);  WriteLn('Y[',i:2,'] =',y[i]);
  END;
  GotoXY(xpos,ypos+n+1);  WriteLn('x     =',x);
  GotoXY(xpos,ypos+n+2);  WriteLn('h     =',h);
  GotoXY(xpos,ypos+n+3);  WriteLn(nmax,' Teilintervalle');
END;

{---------------------------------------------------------------------------}

PROCEDURE Plot (x, y, xmin, xmax, ymin, ymax: REAL; ch: CHAR);

CONST umax=78; vmax=28;                           { Bildschirmbreite/-hoehe }
VAR   u, v: INTEGER;

BEGIN
  u :=      Round((x-xmin)/(xmax-xmin)*umax);
  v := vmax-Round((y-ymin)/(ymax-ymin)*vmax);
  GotoXY(u,v);  Write(ch);
END;

{---------------------------------------------------------------------------}

PROCEDURE RungeKutta (VAR Y, Yo: Vektor; Xstart, Xend: REAL;
                      VAR h: REAL; AbsErr: REAL;
                          n: INTEGER; VAR nmax: INTEGER);

VAR C1, C2, C3, C4, z, dY: Vektor;
    x, xo, yalt, Tol: REAL;
    j, k: INTEGER;

BEGIN
  yalt := 0; nmax := 4;
  REPEAT
    h := ABS((Xend-Xstart)/nmax);                            { Schrittweite }
    FOR k := 1 TO n DO Y[k] := Yo[k];                  { Anfangsbedingungen }
    FOR j := 0 TO Pred(nmax) DO                     { Runge-Kutta-Iteration }
    BEGIN
      xo := Xstart+j*h;
      DGL(dY, Y, xo);
      FOR k := 1 TO n DO
      BEGIN
        C1[k] := h*dY[k];
        z[k] := Y[k]+C1[k]/2;
      END;
      DGL(dY, z, (xo+h/2));
      FOR k := 1 TO n DO
      BEGIN
        C2[k] := h*dY[k];
        z[k] := Y[k]+C2[k]/2;
      END;
      DGL(dY, z, (xo+h/2));
      FOR k := 1 TO n DO
      BEGIN
        C3[k] := h*dY[k];
        z[k] := Y[k]+C3[k];
      END;
      DGL(dY, z, (xo+h));
      FOR k := 1 TO n DO
      BEGIN
        C4[k] := h*dY[k];
        y[k] := y[k]+(C1[k]+2*C2[k]+2*C3[k]+C4[k])/6;
      END;
    END;
    Tol := ABS(yalt-y[1]);                               { Abbruchkriterium }
    nmax := nmax*2;                                         { Intervallzahl }
    yalt := Y[1];
  UNTIL Tol <= AbsErr;
END;

{---------------------------------------------------------------------------}

BEGIN { Partial_DGL_Solver }
                                                      { Anfangsbedingungen: }
  n     := 10;                                     { Zahl der Stuetzstellen }
  t0    := 0.0;                    { Zeit }          { t-Anfangsbedingungen }
  tend  := 1200.0;
  dt    := 125.0;
  Xstart:= 0.0;                    { Laenge }        { x-Anfangsbedingungen }
  Xend  := 100.0;
  dx    := (Xend-Xstart)/(n-1);                            { x-Schrittweite }
  Yo[1] := 0.0;            { Temperatur }            { y-Anfangsbedingungen }
  Yo[n] := 0.0;
  FOR i := 2 TO Pred(n) DO
  BEGIN
    x := Xstart+Pred(i)*dx;                                 { Anfangsprofil }
    Yo[i] := x;
  END;
  AbsErr:= 1.0E-3;                                       { Absoluter Fehler }
  ymin  := 0.0;                                  { y-Ausdehnung der Graphik }
  ymax  := 100.0;

  ClrScr;
  t := t0;
  k := 0;
  WHILE t <= tend DO                         { Zeitlicher Verlauf: t0..tend }
  BEGIN
    RungeKutta(Y, Yo, t0, t, h, AbsErr, n, nmax);       { DGL-System loesen }
    k := Succ(k);
    FOR i := 1 TO n DO                                  { Loesungen plotten }
    BEGIN
      x := Xstart+Pred(i)*dx;
      Plot(x, Y[i], Xstart, Xend, 0, 100, Chr(47+k));    { Zeichen: 0 bis 9 }
    END;
    Ausgabe(Y, x, h, n, nmax);                         { Loesungen ausgeben }
    t := t+dt;
  END;
  GotoXY(0,0);
END.