145 lines
3.2 KiB
ObjectPascal
145 lines
3.2 KiB
ObjectPascal
program GuassianElimination;
|
|
|
|
// Modified from:
|
|
// R. Sureshkumar (10 January 1997)
|
|
// Gregory J. McRae (22 October 1997)
|
|
// http://web.mit.edu/10.001/Web/Course_Notes/Gauss_Pivoting.c
|
|
|
|
{$APPTYPE CONSOLE}
|
|
|
|
{$R *.res}
|
|
|
|
uses
|
|
System.SysUtils;
|
|
|
|
type
|
|
TMatrix = class
|
|
private
|
|
_r, _c : integer;
|
|
data : array of TDoubleArray;
|
|
function getValue(rIndex, cIndex : integer): double;
|
|
procedure setValue(rIndex, cIndex : integer; value: double);
|
|
public
|
|
constructor Create (r, c : integer);
|
|
destructor Destroy; override;
|
|
|
|
property r : integer read _r;
|
|
property c : integer read _c;
|
|
property value[rIndex, cIndex: integer]: double read getValue write setValue; default;
|
|
end;
|
|
|
|
|
|
constructor TMatrix.Create (r, c : integer);
|
|
begin
|
|
inherited Create;
|
|
self.r := r; self.c := c;
|
|
setLength (data, r, c);
|
|
end;
|
|
|
|
destructor TMatrix.Destroy;
|
|
begin
|
|
data := nil;
|
|
inherited;
|
|
end;
|
|
|
|
function TMatrix.getValue(rIndex, cIndex: Integer): double;
|
|
begin
|
|
Result := data[rIndex-1, cIndex-1]; // 1-based array
|
|
end;
|
|
|
|
procedure TMatrix.setValue(rIndex, cIndex : integer; value: double);
|
|
begin
|
|
data[rIndex-1, cIndex-1] := value; // 1-based array
|
|
end;
|
|
|
|
// Solve A x = b
|
|
procedure gauss (A, b, x : TMatrix);
|
|
var rowx : integer;
|
|
i, j, k, n, m : integer;
|
|
amax, xfac, temp, temp1 : double;
|
|
begin
|
|
rowx := 0; // Keep count of the row interchanges
|
|
n := A.r;
|
|
for k := 1 to n - 1 do
|
|
begin
|
|
amax := abs (A[k,k]);
|
|
m := k;
|
|
// Find the row with largest pivot
|
|
for i := k + 1 to n do
|
|
begin
|
|
xfac := abs (A[i,k]);
|
|
if xfac > amax then
|
|
begin
|
|
amax := xfac;
|
|
m := i;
|
|
end;
|
|
end;
|
|
|
|
if m <> k then
|
|
begin // Row interchanges
|
|
rowx := rowx+1;
|
|
temp1 := b[k,1];
|
|
b[k,1] := b[m,1];
|
|
b[m,1] := temp1;
|
|
for j := k to n do
|
|
begin
|
|
temp := a[k,j];
|
|
a[k,j] := a[m,j];
|
|
a[m,j] := temp;
|
|
end;
|
|
end;
|
|
|
|
for i := k+1 to n do
|
|
begin
|
|
xfac := a[i, k]/a[k, k];
|
|
for j := k+1 to n do
|
|
a[i,j] := a[i,j]-xfac*a[k,j];
|
|
b[i,1] := b[i,1] - xfac*b[k,1]
|
|
end;
|
|
end;
|
|
|
|
// Back substitution
|
|
for j := 1 to n do
|
|
begin
|
|
k := n-j + 1;
|
|
x[k,1] := b[k,1];
|
|
for i := k+1 to n do
|
|
begin
|
|
x[k,1] := x[k,1] - a[k,i]*x[i,1];
|
|
end;
|
|
x[k,1] := x[k,1]/a[k,k];
|
|
end;
|
|
end;
|
|
|
|
|
|
var A, b, x : TMatrix;
|
|
|
|
begin
|
|
try
|
|
// Could have been done with simple arrays rather than a specific TMatrix class
|
|
A := TMatrix.Create (4,4);
|
|
// Note ideal but use TMatrix to define the vectors as well
|
|
b := TMatrix.Create (4,1);
|
|
x := TMatrix.Create (4,1);
|
|
|
|
A[1,1] := 2; A[1,2] := 1; A[1,3] := 0; A[1,4] := 0;
|
|
A[2,1] := 1; A[2,2] := 1; A[2,3] := 1; A[2,4] := 0;
|
|
A[3,1] := 0; A[3,2] := 1; A[3,3] := 2; A[3,4] := 1;
|
|
A[4,1] := 0; A[3,2] := 0; A[4,3] := 1; A[4,4] := 2;
|
|
|
|
b[1,1] := 2; b[2,1] := 1; b[3,1] := 4; b[4,1] := 8;
|
|
|
|
gauss (A, b, x);
|
|
|
|
writeln (x[1,1]:5:2);
|
|
writeln (x[2,1]:5:2);
|
|
writeln (x[3,1]:5:2);
|
|
writeln (x[4,1]:5:2);
|
|
|
|
readln;
|
|
except
|
|
on E: Exception do
|
|
Writeln(E.ClassName, ': ', E.Message);
|
|
end;
|
|
end.
|