where : ibrtses delphi

Delphi - solving equations

disclaimer

the source code of this page may not appear correctly in certain browsers
due to special characters. Have a look at the source of this HTML page

Theory

```equations of the type N equations & N unknown can be solved
with the gaussian reduction algorithm.

a11*x1+a12*x2+......+a1n*xn=c1
a21*x1+a22*x2+......+a2n*xn=c2
.
.
.
an1*x1+an2*x2+......+ann*xn=cn

Fill them into a matrix as :

|a11 a12 ..... a1n|   |x1|   |c1|
|a21 a22 ..... a2n|   |x2|   |c2|
|.                | * |. | = |. |
|.                |   |. |   |. |
|.                |   |. |   |. |
|an1 an2 ..... ann|   |xn|   |cn|

where the
aij are the constant coefficients,
xi the unknown,
ci again constants

Whole rows, including aij and ci are interchangeable, leave the x.
The row with highest absolute value of a?1 is taken as top row.
by dividing the row with the leading a, the first aij becomes 1.

Now, from the second row down, take the one with the highest absolute
value and take it as second row. repeat the above steps

Until the matrix looks like :

|1 * ..... *|   |x1|   |*|
|0 1 ..... *|   |x2|   |*|
|0 0 1      | * |. | = |*|
|0 0        |   |. |   |*|
|0 0        |   |. |   |*|
|0 0 0 ... 1|   |xn|   |*|

Now as the solution of Xn is known, the lines can be added
from the bottom up in the same weighted manner until the matrix
looks like :

|1 0 0 0 . 0|   |x1|   |*|
|0 1 0 0 . 0|   |x2|   |*|
|0 0 1 0 . 0| * |. | = |*|
|0 0 0 1 0 0|   |. |   |*|
|0 0 0 0 1 0|   |. |   |*|
|0 0 0 0 0 1|   |xn|   |*|

The constant terms * are the solutions. The top * is x1 and so on.

Now the program :

```

Implementation

``` class TLinEq :
solves linear equations :

|m11[0,0]   m12[0,1]   m13[0,2]  m14[0,3]   .. m1N[0,N-1]   | | x1 |    |c1[0]  |
|m21[1,0]   m22[1,1]   m23[1,2]  m24[1,3]   .. m2N[1,N-1]   | | x2 |    |c2[1]  |
|m31[2,0]   m32[2,1]   m33[2,2]  m34[2,3]   .. m3N[2,N-1]   | | x3 |    |c3[2]  |
|m41[3,0]   m42[3,1]   m43[3,2]  m44[3,3]   .. m4N[3,N-1]   | | x4 |  = |c4[3]  |
|..         ..         ..        ..         .. ..           | | .. |    | ..    |
|mN1[N-1,0] mN2[N-1,1] mN3[0,2]  mN4[N-1,3] .. mNN[N-1,N-1] | | xN |    |cN[N-1]|

the matix is the property M, the constants are property C.
note the indices are 0 based !!
solve returns true, when the solution was possible.

eg :
var eq:TLinEq;

eq:=TLinEq.create(2);
eq.M[0,0]:=2; eq.M[0,1]:=1; eq.C[0]:=3;
eq.M[1,0]:=1; eq.M[1,1]:=-1;eq.C[1]:=0;

if eq.solve then
begin
result[1]:=eq.C[0];
result[2]:=eq.C[1];
end;

type
TEVect=class(TObject)
private
fp:pointer;
fsize:integer;
function getitem(i:integer):extended;
procedure setitem(i:integer;z:extended);
public
constructor create(size:integer);
constructor createnull(size:integer);
destructor destroy;     override;
procedure clear;
procedure mult(z:extended);  //multiply with z
procedure sub(v:TEVect);     // this:=this-v
function copy:TEVect;
procedure submul(z:extended;v:TEVect);     // this:=this-z*v
published
property Q[index:integer]:extended read getitem write setitem; default;
end;

TLinEq = class(TObject)
private
{ Private declarations }
vec:TList;
equs:integer;
function getMitem(index1,index2:integer):extended;
procedure setMitem(index1,index2:integer;z:extended);
function getCitem(index:integer):extended;
procedure setCitem(index:integer;z:extended);
protected
{ Protected declarations }
function greatestpivot(n:integer):integer;
procedure swaplines(a,b:integer);
public
{ Public declarations }
constructor create(dim:integer);
destructor destroy; override;
function solve:boolean;
published
{ Published declarations }
property M[index1,index2:integer]: extended read getMitem write setMitem;
property C[index:integer]: extended read getCitem write setCitem;
end;
......................................................................
type MyVectException = class(Exception);
type extptr=^extended;
.......................................................................
function TEVect.getitem(i:integer):extended;
var p:extptr;
k:integer absolute p;
begin
if (i < 0)or(i >=fsize) then
raise(MyVectException.create('LinEq:index out of vector boundary'));
p:=fp;
k:=k+i*sizeof(extended);
result:=p^;
end;
procedure TEVect.setitem(i:integer;z:extended);
var p:extptr;
k:integer absolute p;
begin
if (i < 0)or(i >=fsize) then
raise(MyVectException.create('LinEq:index out of vector boundary'));
p:=fp;
k:=k+i*sizeof(extended);
p^:=z;
end;
constructor TEVect.create(size:integer);
begin
inherited create;
fsize:=size;
Getmem(fp,fsize*sizeof(extended));
end;
constructor TEVect.createnull(size:integer);
begin
inherited create;
fsize:=size;
Getmem(fp,fsize*sizeof(extended));
FillChar(fp^,fsize*sizeof(extended),0); // clear
end;
destructor TEVect.destroy;
begin
Freemem(fp,fsize*sizeof(extended));
inherited destroy;
end;
procedure TEVect.clear;
begin
FillChar(fp^,fsize*sizeof(extended),0);
end;
procedure TEVect.mult(z:extended);  //multiply with z
var i:integer;
//    p:extptr;
//    k:integer absolute p;
begin
for i:=0 to fsize-1 do Q[i]:=Q[i]*z;
end;
procedure TEVect.sub(v:TEVect);     // this:=this-v
var i:integer;
begin
for i:=0 to fsize-1 do Q[i]:=Q[i]-v[i];
end;
function TEVect.copy:TEVect;
begin
result:=TEVect.create(fsize);
move(fp^,result.fp^,fsize*sizeof(extended));
end;
procedure TEVect.submul(z:extended;v:TEVect);     // this:=this-z*v
var i:integer;
begin
for i:=0 to fsize-1 do Q[i]:=Q[i]-z*v[i];
end;
......................................................................
function TLinEq.getMitem(index1,index2:integer):extended;
begin
result:=TEVect(vec.items[index1])[index2];
end;
procedure TLinEq.setMitem(index1,index2:integer;z:extended);
begin
TEVect(vec.items[index1])[index2]:=z;
end;
function TLinEq.getCitem(index:integer):extended;
begin
result:=TEVect(vec.items[index])[equs];
end;
procedure TLinEq.setCitem(index:integer;z:extended);
begin
TEVect(vec.items[index])[equs]:=z;
end;
function TLinEq.greatestpivot(n:integer):integer;
var i,j:integer;
t,max:extended;
begin
j:=0;max:=0;
for i:=n to equs-1 do begin
t:=abs(M[i,n]);
if t > max then begin
max:=t;
j:=i;
end;
end;
result:=j;
end;
procedure TLinEq.swaplines(a,b:integer);
var v:TEVect;
begin
v:=vec.items[a];
vec.items[a]:=vec.items[b];
vec.items[b]:=v;
end;
constructor TLinEq.create(dim:integer);
var i:integer;
v:TEVect;
begin
inherited create;
vec:=TList.create;
for i:=1 to dim do begin
v:=TEVect.create(dim+1);
v.clear;
end;
equs:=dim;
end;
destructor TLinEq.destroy;
var i:integer;
v:TEVect;
begin
for i:=0 to equs-1 do begin
v:=vec[i];
v.destroy;
end;
vec.destroy;
inherited destroy;
end;
//
Jean Bernard found an error that arised from copy/pasting :
the current line is the one with the greatest pivot.
//
function TLinEq.solve:boolean;
var l1,l2,u:integer;
t:extended;
begin
l1:=0; result:=true;
while (result)and(l1 < equs) do begin
u:=greatestpivot(l1);
t:=M[u,l1];
if (abs(t) < 1e-60) then result:=false
else begin
if (u < > l1)then swaplines(u,l1);
TEVect(vec.items[l1]).mult(1/t);   // [l1,l1]:=1
for l2:=0 to equs-1 do begin
if (l1 < > l2) then begin
t:=M[l2,l1];
if (abs(t) > 1e-60) then
TEVect(vec.Items[l2]).submul(t,TEVect(vec.items[l1]));
end;
end;
inc(l1);
end;
end;
end;

```

findings

```a standard approach is implemented. 2x2 and 3x3 systems can be
solved faster with a simpler approach without loops.
While rather big systems can be solved, it should be noted
that special approaches exist to solve sparse systems and
for those with bands.
```

Feedback is welcome