Saturday, March 3, 2012

Jacobi method pascal / delphi

In numerical linear algebra, the Jacobi method is an algorithm for determining the solutions of a system of linear equations with largest absolute values in each row and column dominated by the diagonal element. Each diagonal element is solved for, and an approximate value plugged in. The process is then iterated until it converges. This algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization. The method is named after German mathematician Carl Gustav Jakob Jacobi.





Given a square system of n linear equations:
A\mathbf x = \mathbf b
where:
A=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}.
Then A can be decomposed into a diagonal component D, and the remainder R:
A=D+R \qquad \text{where} \qquad D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \qquad R = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ a_{21} & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix} ,\qquad \mathbf{x}^{(k+1)} = D^{-1} (\mathbf{b} - R \mathbf{x}^{(k)}).
The element-based formula is thus:
x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i -\sum_{j\ne i}a_{ij}x^{(k)}_j\right),\quad i=1,2,\ldots,n.
Note that the computation of xi(k+1) requires each element in x(k) except itself. Unlike the Gauss–Seidel method, we can't overwrite xi(k) with xi(k+1), as that value will be needed by the rest of the computation. The minimum amount of storage is two vectors of size n.

So here is the code:

program Project1;
{$APPTYPE CONSOLE}
uses
  SysUtils;
 const e = 0.000001;
var n,i,j,k,m,i1,j1,l:integer; a:array[1..5,1..5]of extended;
b:array[1..5,1..5]of extended; tk:array[1..5,1..5]of extended;
q,p,d,r,c,s,sign,ab,eb:real; t:array[1..5,1..5]of extended;
cm: array[1..5,1..5]of extended;

begin
writeln('Input n '); // Matrix size
read(n);
writeln('Input matrix A ');

// Matrix input
  for i:=1 to n do
    for j:=1 to n do
      read(a[i,j]);

    ab:=a[1,2];
    eb:=1;

    for i1:=1 to n do
      for j1:=1 to n do
        if i1<>j1 then
        begin
          if abs(a[i1,j1])>ab then
            begin
              i:=i1;
              j:=j1
            end;
        end;

  for i1:=1 to n do
  for j1:=1 to n do
    if i1=j1 then
      begin
        t[i1,i1]:=1;
        tk[i1,i1]:=1;
      end
    else
      t[i1,j1]:=0;

while eb>e do
begin           //---1
    p:=2*a[i,j];
    q:=a[i,i]-a[j,j];
    d:=sqrt(p*p+q*q);
    sign:=p*q;

  if sign<>0 then
  sign:=sign/abs(sign)
    else
    sign:=0;
                  //---2
  if (q<>0)and((abs(p))>=(abs(q))) then
    begin
      r:=abs(q)/(2*d);
      c:=sqrt(0.5+r);
      s:=(sqrt(0.5-r))*sign
    end
  else
  if (q<>0)and ((abs(p))<(abs(q))) then
    begin
      r:=abs(q)/(2*d);
      c:=sqrt(0.5+r);
      s:=(abs(p))*sign/(2*c*d)
    end
  else
  begin
    c:=sqrt(2)/2;
    s:=c
  end;
                  // ----3
  b[i,i]:=c*c*a[i,i]+s*s*a[j,j]+2*c*s*a[i,j];
  b[j,j]:=s*s*a[i,i]+c*c*a[j,j]-2*c*s*a[i,j];
  b[i,j]:=0;    //---4
  b[j,i]:=0;

  //----5
  for m:=1 to n do
    if (m<>i)and(m<>j) then
    begin
      b[i,m]:=c*a[m,i]+s*a[m,j];
      b[m,i]:=b[i,m];
      b[j,m]:=-s*a[m,i]+c*a[m,j];
      b[m,j]:=b[j,m];
    end;
  //--6
  for m:=1 to n do
  for l:=1 to n do
    if (m<>i)and(m<>j)and(l<>i)and(l<>j) then
    begin
      b[m,l]:=a[m,l];
      b[l,m]:=b[m,l]
    end;
   //---matrix T
  for i1:=1 to n do
  for j1:=1 to n do
    if i1=j1 then
    begin
      tk[i1,i1]:=1;
      cm[i1,i1]:=0;
    end
    else
  begin
    cm[i1,j1]:=0;
    tk[i1,j1]:=0
  end;

tk[i,i]:=c;
tk[j,j]:=c;
tk[j,i]:=-s;
tk[i,j]:=s;

for i1:=1 to n do
  for j1:=1 to n do
  for k:=1 to n do
    cm[i1,j1]:=cm[i1,j1]+tk[i1,k]*t[k,j1];

for i1:=1 to n do
  for j1:=1 to n do
    begin
      t[i1,j1]:=cm[i1,j1];
    end ;

    eb:=b[1,2];

for i1:=1 to n do
  for j1:=1 to n do
  if i1<>j1 then
    begin
      if abs(b[i1,j1])>eb then
      begin
        eb:=abs(b[i1,j1]);
        i:=i1;
        j:=j1;
      end;
    end;

  eb:=abs(eb);

    for i1:=1 to n do
      for j1:=1 to n do
      a[i1,j1]:=b[i1,j1];

end;
    
writeln('Vlasni Paru');
for i1:=1 to n do
 begin
    for j1:=1 to n do
    write(b[i1,j1]:4:4,' ');
    writeln;
  end;


writeln;
writeln('Vectors');
  for i1:=1 to n do
  begin
    for j1:=1 to n do
    writeln('X',j1,'=', t[i1,j1]:2:4);
    writeln('----------');
  end;
    readln;
    readln;

end.

Posted by at
Labels:

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)