previous overview
program MultiGrid;
{
The End of MultiGrid
====================
but only for 1-D ...
Let the system of equations be given by S.w = b .
Iterate: M := I - T.S ; b := (I + M).b ; S := I - M^2 .
The background of this being (I - M)^(-1) = (I + M)/(I - M^2).
The matrix T is a preconditioner, producing normed equations.
For tri-diagonal 1-D systems, the matrix (I - M^2) recursively
reduces to blocks around the main diagonal. It finally becomes
a diagonal matrix which is normed, giving the Exact solution.
CopyLefted by: Han de Bruijn (HdB)
}
const
NN = 11 ; { Number of unknowns }
var
e : array[1..2,1..2] of double;
s,t : array[1..NN,0..2] of double;
b,p : array[1..NN] of double;
nr,no : array[1..NN] of byte;
effe : text; { LOG file }
L : byte ; getal : double ;
procedure Element(eps : double);
{
Define Finite Element matrix
---------------------------- }
var
a,b : double;
begin
{
With Safety Condition:
a.b = 0.5*(1 + eps)*0.5*(1 - eps)
= (1 - sqr(eps))/4 < 1/4
Because: eps = getal = Random < 1
}
a := 0.5*(1 - eps);
b := 0.5*(1 + eps);
e[1,1] := +a ; e[1,2] := -a;
e[2,1] := -b ; e[2,2] := +b;
end;
procedure Boekhouden(eerst : boolean);
{
Administration
}
var
i : byte;
begin
if eerst then
for i := 1 to NN do
nr[i] := i;
{ Inverse: }
for i := 1 to NN do
no[nr[i]] := i;
end;
procedure Normeren;
{
Make diagonal = unity
}
var
i : byte;
d : double;
begin
for i := 1 to NN do
begin
d := 1/s[i,1];
s[i,0] := s[i,0] * d;
s[i,2] := s[i,2] * d;
s[i,1] := 1;
b[i] := b[i] * d;
end;
end;
procedure Schoonmaken;
{
Clear global matrix and vector
}
var
i : byte;
begin
for i := 1 to NN do
begin
s[i,0] := 0 ;
s[i,1] := 0 ;
s[i,2] := 0 ;
b[i] := 0;
end;
end;
procedure A_Symmetrisch;
{
Fill global matrix & vector
}
var
n,i,j,ii,jj : byte;
begin
for n := 1 to NN-1 do
begin
for i := 1 to 2 do
for j := 1 to 2 do
begin
ii := n+i-1 ; jj := j-i+1 ;
s[ii,jj] := s[ii,jj] + e[i,j];
end;
end;
b[1] := 1;
s[1,1] := 1 ; s[1,2] := 0 ;
b[NN] := 0;
s[NN,0] := 0 ; s[NN,1] := 1 ;
end;
procedure Oplossen;
{
Decomposition of tri-diagonal System
------------------------------------ }
var
k : integer;
diag, pivot : extended;
begin
for k := 1 to NN do
begin
diag := s[k,1];
if diag = 0 then begin
Writeln('Oplossen: 0 on diagonal at: ',k);
Halt;
end ;
pivot := s[k+1,0]/diag ;
s[k+1,0] := pivot;
if pivot = 0 then Continue ;
s[k+1,1] := s[k+1,1] - pivot * s[k,2] ;
end;
end;
procedure Oprollen;
{
Solution of tri-diagonal System
------------------------------- }
var
k : integer;
diag, pivot : extended;
begin
for k := 1 to NN do
begin
pivot := s[k+1,0] ;
if pivot = 0 then Continue ;
b[k+1] := b[k+1] - pivot * b[k] ;
end;
for k := NN downto 1 do
begin
pivot := b[k];
diag := s[k,1];
if diag = 0 then begin
Writeln('Oplossen: 0 on diagonal at: ',k);
Halt;
end ;
pivot := pivot - s[k,2] * b[k+1] ;
b[k] := pivot/diag;
end;
end;
procedure Vullen(term : double);
{
Define System of Equations
-------------------------- }
var
k : integer;
begin
{ Bulk Assembly: }
for k := 2 to NN-1 do
begin
s[k,1] := 1 ; b[k] := 0;
s[k,0] := - 0.5*(1 + term);
s[k,2] := - 0.5*(1 - term);
end;
{ Left boundary: }
s[1,1] := 1;
s[1,2] := 0;
b[1] := 1;
{ Right boundary: }
s[NN,1] := 1;
s[NN,0] := 0;
b[NN] := 0;
end;
procedure Afdrukken;
{
Print out Solution
}
var
i : byte;
begin
for i := 1 to NN do
Write(b[no[i]]:7:4);
Writeln;
end;
procedure Bekijken(onder,boven : byte);
{
Take a snapshot in the LOG
}
var
i : byte;
begin
for i := onder+1 to boven do
Writeln(effe,s[i,0]:9:5,
' ',s[i,1]:9:5,
' ',s[i,2]:9:5);
for i := onder+1 to boven do
Write(effe,b[i]:9:5,' ');
Writeln(effe) ; Writeln(effe);
end;
procedure Newton(onder,boven : byte);
{
Newton-Raphson Multigrid
------------------------ }
var
i,ii : byte;
midden : byte;
d : double;
begin
{
Renormalization
}
for i := onder+1 to boven do
begin
d := 1/s[i,1];
s[i,0] := s[i,0] * d;
s[i,2] := s[i,2] * d;
s[i,1] := 1;
b[i] := b[i] * d;
end;
Bekijken(onder,boven);
if (onder+1 = boven) then Exit; { We are done ! }
{
M := I - S ; b := (I + M) * b
}
p := b;
for i := onder+1 to boven do
begin
p[i] := b[i] - s[i,0]*b[i-1] - s[i,2]*b[i+1];
end;
b := p;
{
(I + M).(I - M) = I - M^2 -> main diagonal
}
for i := onder+1 to boven do
s[i,1] := 1 - s[i,0]*s[i-1,2] - s[i,2]*s[i+1,0];
{
(I + M).(I - M) = I - M^2 -> off diagonal terms
}
{ Two cases: # unknowns even or odd }
if ((boven - onder) mod 2) = 0 then
midden := onder + ((boven - onder) div 2);
if ((boven - onder) mod 2) = 1 then
midden := onder + ((boven - onder + 1) div 2);
{ Block out Odd indices: }
for i := onder+1 to midden do
begin
ii := onder + 2*(i-onder) - 1;
t[i,1] := s[ii,1];
t[i,0] := - s[ii,0]*s[ii-1,0];
t[i,2] := - s[ii,2]*s[ii+1,2];
p[i] := b[ii];
no[i] := nr[ii];
end;
t[onder+1,0] := 0;
t[midden,2] := 0;
{ Block out Even indices: }
for i := midden+1 to boven do
begin
ii := onder + 2*(i-midden);
t[i,1] := s[ii,1];
t[i,0] := - s[ii,0]*s[ii-1,0];
t[i,2] := - s[ii,2]*s[ii+1,2];
p[i] := b[ii];
no[i] := nr[ii];
end;
t[midden+1,0] := 0;
t[boven,2] := 0;
{ Recursively: }
s := t ; b := p ;
nr := no ;
Newton(onder,midden);
Newton(midden,boven);
end;
procedure FDM(eps : double);
{
Finite Difference Method
------------------------ }
begin
{ Conventional }
Boekhouden(true);
Vullen(eps);
Oplossen;
Oprollen;
Afdrukken;
{ Newton-Raphson }
Vullen(eps);
Newton(0,NN);
Boekhouden(false);
Afdrukken;
end;
procedure FEM(eps : double);
{
Finite Element Method
--------------------- }
begin
{ Conventional }
Boekhouden(true);
Element(eps);
Schoonmaken;
A_Symmetrisch;
Oplossen;
Oprollen;
Afdrukken;
{ Newton-Raphson }
Schoonmaken;
A_Symmetrisch;
Newton(0,NN);
Boekhouden(false);
Afdrukken;
end;
begin
{
The method should work for any a-symmetric,
tri-diagonal and "safe" system of equations
}
Assign(effe,'newton.log');
Rewrite(effe);
for L := 1 to 4 do
begin
Writeln;
getal := Random;
FDM(getal);
FEM(getal);
end;
Close(effe);
end.