alg072.ma

数值计算方法的Mathmatic实现

(* GAUSS-SEIDEL ITERATIVE TECHNIQUE ALGORITHM 7.2
*
*  To solve Ax = b given an initial approximation x(0)
*
*   Input:  the number of equations and unknowns n;
*            the entries A(i,J), 1<=i, J<=n of the matrix A;
*	     The entries B(i), 1<=i<=n of the inhomogeneous term b;
*	     the entries XO(i), 1<=i<=n, of x(0);
*	     tolerance TOL; maximum number of iterations NN.
*
*  Output:  the approximate solution X(1), ..., X(n) or a message
*	     that the maximum number of iterations was exceeded.
*)
Print["A(1,1),A(1,2)..A(1,n+1)\n"];
Print["A(2,1),A(2,2)..A(2,n+1)\n"];
Print["A(n,1),A(n,2)..A(n,n+1)\n"];
Print["initial approximation\n"];
Print["Separate entries with blank\n"];
AA = InputString["This is Gauss-Seidel Method for Linear Systems.\n
The array will be input from a text file in the order:(see screen)\n
Has the input file been created?  Enter 'yes' or 'no'\n"];
OK = 0;
If[AA == "yes" || AA == "y" || AA == "Y",
   NAME=InputString["Input the file name in the form - \n
      drive:\ name.ext      for example\n
         A:\\DATA.DTA\n"];
   INP = OpenRead[NAME];
   OK = 0;
   While[OK == 0,
      n = Input["Input the number of equations n - an integer\n"];
      If[n > 0,
	 For[i = 1,i <= n,i++,
	    For[J = 1,J <= n+1,J++,
	       A[i-1,J-1] = Read[INP,Number];
   	    ];
	 ];
(* Use X1 for X0 *)
	 For[i = 1, i <= n, i++,
	    X1[i-1] = Read[INP,Number];
	 ];
         OK = 1;
	 Close[INP],
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      TOL = Input["Input the tolerance\n"];
      If[TOL <= 0,
         Input["Tolerance must be positive.\n
         Enter 1 [enter] to continue\n"],
         OK = 1;
      ];
   ];
   OK = 0;
   While[OK == 0,
      NN = Input["Input maximum number of iterations\n"];
      If[NN <= 0,
         Input["Must be a positive integer.\n
         Enter 1 [enter] to continue\n"],
         OK = 1
      ];
   ],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   (* Step 1 *)
   K = 1;
   OK = 0;
   (* Step 2 *)
   While[OK == 0 && K <= NN,
      (* ERR is used to test accuracy - it measures the infinity-norm *)
      ERR=0;
      (* Step 3 *)
      For[i = 1, i <= n, i++,
	 S = 0;
	 For[J = 1, J <= n, J++,
	    S = S-A[i-1,J-1]*X1[J-1];
	 ];
	 S = (S+A[i-1,n])/A[i-1,i-1];
	 If[Abs[S] > ERR,
	    ERR = Abs[S];
	 ];
	 X1[i-1] = X1[i-1]+S;
      ];
      (* Step 4 *)
      If[ERR <= TOL,
      (* Step 5 - Process is complete *)
	 OK = 1,
      (* Step 6 - is not used since only one vector is required *)
      K = K+1;
     ];
   ];
   If[OK == 0,
      Write[OUP,"Maximum number of iterations exceeded\n"],
      (* Step 7 - Procedure completed unsuccessfully *)
      FLAG = Input["Select output destination\n
        1. Screen\n
        2. Text file\n
        Enter 1 or 2\n"];
      If[FLAG == 2,
         NAME = InputString["Input the file name\n
            For example:   output.dta\n"];
         OUP = OpenWrite[NAME,FormatType->OutputForm],
         OUP = "stdout";
      ];
      Write[OUP,"GAUSS-SEIDEL METHOD FOR LINEAR SYSTEMS\n"];
      Write[OUP,"\n"];
      Write[OUP,"The solution vector is:\n"];
      For[i = 1, i <= n, i++,
         Write[OUP,N[X1[i-1],9]];
      ];
      Write[OUP,"\n"];
      Write[OUP,"Using ",K," iterations\n"];
      Write[OUP,"with tolerance ",TOL," in infinity-norm\n"];
      If[OUP == "OutputStream[",NAME," 3]",
         Print["Output file: ",NAME," created successfully\n"];
         Close[OUP];
      ];
   ];
];