Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Object Pascal: кроссплатформенные технологии > Метод Якоби

Автор: Mr.Picwick 11.12.2004, 00:30

Привет всем, дикость конечно, но мне надо сделать лабу на Паскале.
Тема: Решение СЛАУ методом Якоби. Прогу я написал, но выдается ERROR
Паскаль 7.1 Win2K, Celeron433.

Всем спасибо.

Код


  program ya;
type
  vector = array[1..4] of real;
  matrix = array[1..4,1..4] of real;
const
  b : vector = (6,-12,1,3);
  a : matrix = ((3,1,-1,2),(-5,1,3,-4),(2,0,1,-1),(1,-5,3,-3));
  n : integer = 4;
  eps : real = 1;
var
  i,j,max,k : integer;
  sum1,sum2 : real;
  beta,x_curr,x_prev,x_tmp  : vector;
  alpha : matrix;
begin

  {alpha calculating}
  {for i:=1 to n do
   for j:=1 to n do
     alpha[i,j]:=-a[i,j]/a[i,i];}

  {beta calculating}
  for i:=1 to n do
    beta[i]:=b[i]/a[i,i];

  {initial x calc}
  for i:=1 to n do
   x_prev[i]:=beta[i];

  sum1:=0;
  sum2:=0;
  k:=0;
  repeat
    inc(k);
    {x_tmp to x_prev}
     if k <> 1 then
       for i:=1 to n do
         x_prev[i]:=x_tmp[i];

    {x_current  calculating}
    for i:=1 to n do
      begin
        sum1:=0;
        {first sum calc}
        for j:=1 to i-1 do
         sum1:=sum1 + a[i,j] * x_prev[j];

   {!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Здесь выдается ошибка Floating-point error}



        sum2:=0;
        {second sum calc}
        for j:=i+1 to n do
         sum2:=sum2 + a[i,j] * x_prev[j];
        {current xi calc}
        x_curr[i]:=(b[i] - sum1 - sum2)/a[i,i];
      end;

    {x_curr to x_tmp}
    for i:=1 to n do
      x_tmp[i]:=x_curr[i];

    {max (x_curr-x_prev) calc}
    max:=1;
    for i:=1 to n-1 do
      if abs(x_curr[max]-x_prev[max]) < abs(x_curr[i+1]-x_prev[i+1]) then
        max:=i+1;

  until  abs(x_curr[max]-x_prev[max]) < eps;

  for i:=1 to n do
   writeln('x[',i,']=',x_curr[i]);
end.

Автор: Гость_nero_schwarz 11.12.2004, 02:44

Цитата

for j:=1 to i-1 do
sum1:=sum1 + a[i,j] * x_prev[j];

я конечно не спец, но в первом цикле, когда у тебя i=1, получается, что j должен изменятся от 1 до 0! и члена a[0,1] массивa у тебя тоже не существует!!!

Автор: Mr.Picwick 11.12.2004, 09:26

все правильно, прога просто не заходит в этот цикл первый раз, остальные разы заходит, это и в алгоритме написано, чтобы первый раз
в цикл суммы не входить.

Может ЭТО зависеть от проца, винды, ТП ?

Автор: Ruslan_ 11.12.2004, 19:36

Я давно писал эту программку. Но помню, что работала исправно. Это точно!

Код


uses crt;
const
  e=0.001;
  num1=5;
  num2=3;

type
  TX=array [1..100] of real;
  TA=array [1..100] of TX;
var a:TA;
    x,tmp:TX;
    n,i:integer;
    infofil:text;

    t:real;
{===================== Checking for file exist ===========================}
function CheckFileExist:boolean;
var z:char;
Begin
    {$I-}Reset(input);{$I+}
    If IoResult<>0 then CheckFileExist:=false
    else CheckFileExist:=true
End;

procedure LoadFromFile(var a:TA;var n:integer);
var i,j:integer;
Begin
  i:=1; n:=-1;
  while not eof(input)and (i<31) do
  begin
    j:=1;
    While not eoln(input)and not eof(input) and (j<32) do
    begin
      Read(a[i][j]);
      inc(j)
    end;
    Readln;
    inc(i);
  end;
  n:=i-1;
End;

Procedure SaveToFile(var x:TX; n:integer);
var i,j:integer;
Begin
  For i:=1 to n do
    Writeln(x[i]:num1:num2);
end;

function CheckMatrix(var a: TA; n:integer):boolean;
var i,j:integer; Aij:real;
begin
  for i:=1 to n do
  begin
    Aij:=0;
    for j:=1 to n do
      if (i<>j) then Aij:=Aij+a[i,j];
    if abs(a[i,i])<Aij then begin CheckMatrix:=false; EXIT; end;
  end;
  CheckMatrix:=true;
end;

procedure Method_Yakobi(var a:TA;var x:TX;n:integer;const e:real);
var S1,S2:real; i,j:integer; xnxt,dx:TX;

function eps:boolean;
var i:integer;
begin
  for i:=1 to n do
  begin
    if abs(xnxt[i]-x[i])<=e then eps:=true else eps:=false
  end;
end;

begin
  FillChar(x,n,'0'); i:=-1;
repeat
  if i>0 then x:=xnxt;
  for i:=1 to n do
  begin
    S1:=0; S2:=0;
    for j:=1 to n do
    begin
      if j<i then S1:=S1+a[i,j]*x[j];
      if j>i then S2:=S2+a[i,j]*x[j]
    end;
    xnxt[i]:=( a[i,n+1]-S1-S2 ) / a[i,i]
  end;
until eps;
x:=xnxt
end;

procedure MultiMatrix(var A: TA; var B,C: TX; n:integer);
var i,j:integer; s:real;
begin
  for i:=1 to n do
  begin
    s:=0;
    for j:=1 to n do
      s:=s+A[i,j]*B[j];
    C[i]:=s;
  end
end;

Procedure SpecSubtractionMatrix(var A:TX;var B: TA;var C:TX; n:integer);
var i:integer;
begin
  for i:=1 to n do
    C[i]:=A[i]-B[i,n+1];
end;


BEGIN
  TextBackGround(0);TextColor(7);ClrScr;

  Assign(input,'Input.txt'); Assign(output,'Output.txt');
  Assign(infofil,'Info.txt');
  Rewrite(infofil);  Rewrite(output);
  Writeln(infofil,'Program was started');
  if CheckFileExist then
  begin
    Reset(input);
    LoadFromFile(a,n)
  end
  else
  begin
    Writeln(infofil,'Error: File Input.txt not found.');
    Close(infofil);
    Close(output);
    EXIT
  end;
  if not CheckMatrix(a,n) then
  begin
    Writeln(infofil,'Error: Bad Matrix.');
    Close(infofil);
    Close(output);
    EXIT;
  end;
  Method_Yakobi(a,x,n,e);
  SaveToFile(x,n);
  MultiMatrix(a,x,tmp,n);
  SpecSubtractionMatrix(tmp,a,tmp,n);
  Writeln; Writeln('Wcheck-up (~f-f):');
  SaveToFile(tmp,n);
  Close(output);
  Close(input);
  Writeln(infofil,'Program successful complete');
  Close(infofil);
END.

Работает с файлом Input.txt и Output.txt. Данные вводятся перечислением их в файле. Строчка - данные для i-го у-ия. В конце каждой строчки, по-моему, должен стоять пробел (не обязательно).
Добавлено @ 19:43
Твоя прога не работает - ошибка 205: Floating point overflow.
Извини, но разбираться времени нет.

Автор: Zero 11.12.2004, 21:07

Цитата(Mr @ 11.12.2004, 00:30)

ошибка Floating-point error

Нас сколько мне помнится это переполнение, для примера можеш попробывать перед циклом напиши

Код


...
sum1:=int(sum1)
...

Ошибки не будет...

Автор: Mr.Picwick 12.12.2004, 14:40

Спасибо чуваки, я знал что я не один!

Проблему я решил - решение чисто математическое, не все СЛАУ можно решить методом Якоби, нужно чтобы система удовлетворяла условиям сходимости, а мой пример НЕ СХОДИТСЯ. Другой пример(СЛАУ) решается очень точно.

Всем спасибо еще раз.