home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / pascal / library / dos / nrpas / tridag.dem < prev    next >
Text File  |  1994-04-11  |  2KB  |  62 lines
  1. PROGRAM d2r4(input,output,dfile);
  2. (* driver for routine TRIDAG *)
  3. LABEL 10,99;
  4. CONST
  5.    np=20;
  6. TYPE
  7.    glnarray = ARRAY [1..np] OF real;
  8. VAR
  9.    k,n: integer;
  10.    diag,superd,subd,rhs,u: glnarray;
  11.    dfile : text;
  12.  
  13. (*$I MODFILE.PAS *)
  14. (*$I TRIDAG.PAS *)
  15.  
  16. BEGIN
  17.    glopen(dfile,'matrx2.dat');
  18. 10:   readln(dfile);
  19.    readln(dfile);
  20.    readln(dfile,n);
  21.    readln(dfile);
  22.    FOR k := 1 to n-1 DO read(dfile,diag[k]);
  23.    readln(dfile,diag[n]);
  24.    readln(dfile);
  25.    FOR k := 1 to n-2 DO read(dfile,superd[k]);
  26.    readln(dfile,superd[n-1]);
  27.    readln(dfile);
  28.    FOR k := 2 to n-1 DO read(dfile,subd[k]);
  29.    readln(dfile,subd[n]);
  30.    readln(dfile);
  31.    FOR k := 1 to n-1 DO read(dfile,rhs[k]);
  32.    readln(dfile,rhs[n]);
  33. (* carry out solution *)
  34.    tridag(subd,diag,superd,rhs,u,n);
  35.    writeln ('the solution vector is:');
  36.    FOR k := 1 to n-1 DO write(u[k]:12:6);
  37.    writeln(u[n]:12:6);
  38. (* test solution *)
  39.    writeln ('(matrix)*(sol''n vector) should be:');
  40.    FOR k := 1 to n-1 DO write(rhs[k]:12:6);
  41.    writeln(rhs[n]:12:6);
  42.    writeln ('actual result is:');
  43.    FOR k := 1 to n DO BEGIN
  44.       IF (k = 1) THEN BEGIN
  45.          rhs[k] := diag[1]*u[1] + superd[1]*u[2]
  46.       END ELSE IF (k = n) THEN BEGIN
  47.          rhs[k] := subd[n]*u[n-1] + diag[n]*u[n]
  48.       END ELSE BEGIN
  49.          rhs[k] := subd[k]*u[k-1] + diag[k]*u[k]
  50.          + superd[k]*u[k+1]
  51.       END
  52.    END;
  53.    FOR k := 1 to n-1 DO write(rhs[k]:12:6);
  54.    writeln(rhs[n]:12:6);
  55.    writeln ('***********************************');
  56.    IF eof(dfile) THEN GOTO 99;
  57.    writeln ('press return for next problem:');
  58.    readln;
  59.    GOTO 10;
  60. 99:   close(dfile)
  61. END.
  62.