Re: 4x4 matrices using Cramer's Rule

Uinseann wrote:

Anyone got any ideas on how to solve 4x4 matrices by using Cramer's
Rule. I've looked in a multitude of math books and surfed the web for
hours to no avail. I can do a 3x3 no problem but unfortunately I dont
seem to be able to see how to do a 4X4. Any advice that anyone might
be able to offer regarding the problem below would be greatly
appreciated.

13 10 0 0 i1 6
-10 13 0 -3 X i2 = 10
0 0 18 -3 i3 0
0 -3 -3 6 i4 5

The following is in Delphi Pascal. I hope it's so much readable that you
can translate it into you own favorite programming language. The outcome
is, iff I've made no mistakes:

-2.38532110091743E-0001 = i1
9.10091743119266E-0001 = i2
2.34250764525994E-0001 = i3
1.40550458715596E+0000 = i4

Han de Bruijn
=====================================================================
Unit Algemeen;
{
General purpose routines
========================
}
INTERFACE

type
integers = array of integer;
matrix = array of array of double;

function Faculteit(n : integer) : integer; { Factorial }
{ Permutations in lexographic order. Algorithm by Edsger W. Dijkstra }
function Dijkstra(waarde : integers; var teken : integer) : integers;
procedure Edsger(N : integer); { Just a Test }

type
VektorRuimte = class(TObject) { Vector Space }
private
NDM : integer; { Dimension }
public
procedure Dimensie(grootte : integer); { Define Dimension }
function Determinant(getallen : matrix) : double;
function DirekteInverse(var getallen : matrix) : boolean;
end;

procedure DirectDemo(NDM : integer);
procedure Uinseann;

IMPLEMENTATION

function als(prop : boolean) : integer;
{ ===
IF like in APL
-------------- }
begin
als := 0;
if prop then als := 1;
end;

function Faculteit;
{
Factorial n!
------------ }
var
k,m : integer;
begin
m := 1;
for k := 2 to n do
m := m*k;
Faculteit := m;
end;

function Dijkstra;
{
Should produce a Permutation (waarde)
in alphabetical (lexicographic) order.
Start with the identical permutation.

According to: E.W.Dijkstra,
"A discipline of Programming",
Prentice Hall, 1976
}
var
i,j,N : integer;

procedure Wissel(k,L : integer);
var
h : integer;
begin
h := waarde[L];
waarde[L] := waarde[k];
waarde[k] := h;
teken := - teken;
end;

begin
N := Length(waarde);

i := N-1;
while (waarde[i-1] >= waarde[i]) do i := i-1;
j := N;
while (waarde[j-1] <= waarde[i-1]) do j := j-1;

Wissel(i-1,j-1);
i := i+1; j := N;

while (i < j) do
begin
Wissel(i-1,j-1);
i := i+1;
j := j-1;
end;

Dijkstra := waarde;
end;

procedure Edsger;
{
Test Permutations
}
var
orde,F,k,i,teken : integer;
rij : integers;

function Extra : string;
var
s : string;
begin
s := '+ ';
if teken < 0 then s := '- ';
Extra := s;
end;

begin
orde := N;
SetLength(rij,orde);
F := Faculteit(orde);
for i := 0 to orde-1 do
rij[i] := i+1;
teken := +1;
for k := 1 to F do
begin
if k > 1 then
rij := Dijkstra(rij,teken);
Write(Extra);
for i := 0 to orde-1 do
Write(rij[i]:3);
Writeln;
end;
end;

function subdet(ii,jj,ndm : integer; getallen : matrix) : double;
{ ======
Subdeterminant of small matrix
------------------------------ }
var
rij : integers;
fout : boolean;
det, term : double;
m,i,j,k, even : integer;

function teken(L : integer) : integer;
begin
teken := 1-2*(L mod 2);
end;

begin
subdet := 0;
fout := (ndm < 1) or (ndm > 8);
if fout then
begin
Writeln('subdet: input Error');
Exit;
end;

SetLength(getallen,ndm+1,ndm+1);
SetLength(rij,ndm-1);
for k := 0 to ndm-2 do
rij[k] := k+1;

det := 0;
even := +1;
for m := 1 to Faculteit(ndm-1) do
begin
if m > 1 then
rij := Dijkstra(rij,even);
term := even;
for k := 1 to ndm-1 do
begin
i := k + als(k >= ii);
j := rij[k-1] + als(rij[k-1] >= jj);
term := term*getallen[i,j];
end;
det := det + term;
end;
subdet := teken(ii+jj) * det;
end;

function VektorRuimte.Determinant(getallen : matrix) : double;
{ ===========
Determinant of small matrix
--------------------------- }
var
rij : integers;
fout : boolean;
det, term : double;
m,i,j,k, even : integer;

function teken(L : integer) : integer;
begin
teken := 1-2*(L mod 2);
end;

begin
Determinant := 0;
fout := (NDM < 1) or (NDM > 8);
if fout then
begin
Writeln('Determinant: input Error');
Exit;
end;

SetLength(getallen,NDM+1,NDM+1);
SetLength(rij,NDM);
for k := 0 to NDM-1 do
rij[k] := k+1;

det := 0;
even := +1;
for m := 1 to Faculteit(NDM) do
begin
if m > 1 then
rij := Dijkstra(rij,even);
term := even;
for k := 1 to NDM do
begin
i := k; j := rij[k-1];
term := term*getallen[i,j];
end;
det := det + term;
end;
Determinant := det;
end;

procedure VektorRuimte.Dimensie;
{
Dimension of Vector Space
}
begin
NDM := grootte;
end;

function VektorRuimte.DirekteInverse;
{
Direct inversion of small matrices
---------------------------------- }
var
i,j : integer;
invert : matrix;
det : double;
begin
DirekteInverse := false;
if ndm < 1 then
begin
Writeln('DirekteInverse: wrong input');
Exit;
end;
SetLength(invert,ndm+1,ndm+1);

for i := 1 to ndm do
for j := 1 to ndm do
invert[i,j] := subdet(j,i,ndm,getallen);

det := 0;
for i := 1 to ndm do
det := det + getallen[i,1]*invert[1,i];
if det = 0 then
begin
Writeln('DirekteInverse: singular matrix');
Exit;
end;

for i := 1 to ndm do
for j := 1 to ndm do
invert[i,j] := invert[i,j]/det;

getallen := invert;
DirekteInverse := true;
end;

procedure Scherm(ndm : integer; beeld : matrix);
{
Print
}
var
k,L : integer;
begin
for L := 1 to ndm do
begin
for k := 1 to ndm do
Write(' ',beeld[k][L]:8:5);
Writeln;
end;
Writeln;
end;

procedure DirectDemo;
{ ==========
Direct inversion
----------------
demonstration program
--------------------- }
var
a,b,c : matrix;
i,j,k : integer;
h : double;
ruimte : VektorRuimte;
begin
SetLength(a,NDM+1,NDM+1);
SetLength(b,NDM+1,NDM+1);
SetLength(c,NDM+1,NDM+1);

{ Random matrix }
for i := 1 to NDM do
for j := 1 to NDM do
a[i,j] := Random;
Scherm(NDM,a);
b := a;

ruimte := VektorRuimte.Create;
ruimte.Dimensie(NDM);
{ Direct inversion }
ruimte.DirekteInverse(b);
Scherm(NDM,b);

{ Multiply original with inverse }
for i := 1 to NDM do
begin
for j := 1 to NDM do
begin
h := 0.;
for k := 1 to NDM do
h := h + a[i,k]*b[k,j];
c[i,j] := h;
end;
end;

{ Result must be unity }
Scherm(NDM,c);
end;

procedure Uinseann;
{ ========
Direct inversion
----------------
demonstration program
--------------------- }
const
NDM : integer = 4;
b : array[1..4] of double = (6,10,0,5);
var
a : matrix;
c : array[0..4] of double;
i,j : integer;
h : double;
ruimte : VektorRuimte;
begin
SetLength(a,NDM+1,NDM+1);
a[1,1] := 13; a[1,2] := 10; a[1,3] := 0; a[1,4] := 0;
a[2,1] := -10; a[2,2] := 13; a[2,3] := 0; a[2,4] := -3;
a[3,1] := 0; a[3,2] := 0; a[3,3] := 18; a[3,4] := -3;
a[4,1] := 0; a[4,2] := -3; a[4,3] := -3; a[4,4] := 6;
Scherm(NDM,a);

ruimte := VektorRuimte.Create;
ruimte.Dimensie(NDM);
{ Direct inversion }
ruimte.DirekteInverse(a);
Scherm(NDM,a);
{ Multiply inverse with load }
for i := 1 to NDM do
begin
h := 0.;
for j := 1 to NDM do
begin
h := h + a[i,j]*b[j];
end;
c[i] := h;
end;

for i := 1 to 4 do
Writeln(c[i]);
end;
{
program proef;

Uses Algemeen;

begin
Uinseann;
end.
}
END.

.