Решение СЛАУ методоми простых итераций и Якоби

chebog · 20.05.2008

Люди, нужно сессию закрыть, всю ночь сегодня не спал. Голова звенит...
В общем, делаю зачетное задание (как всегда в последние минуты...)

Написал код, который по сути должен решать СЛАУ методами простых итераций и Якоби, но когда пытаюсь решить уравнения, ясно видно, что методы не сходятся... Весь вопрос - почему. По всем условиям должны сходиться. Сейчас отлаживаю, но сдавать уже через 40 минут, так что хотелось бы подключить общественность.

Вот код:

Файл iteracii.cpp:
[codebox]// iteracii.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "Classes.h"
#include <math.h>

using namespace std;

void SIM(int size) //Метод простых итераций
{
int i = 0;
double eps;

Vector xip(size);
Vector xi(size);
Vector b(size);
Matrix a(size);
Matrix e(size);

double tau=1; //xip = (a + tau*e)*xi + tau*b;
//последовательность сходится, если каждий из диагональных элементов матрицы по модулю больше суммы модулей всех
//недиагональных элементов соответствующего столбца(строки).

cout<<"Метод простых итераций\n";
e.FillE();
cout<<"Введите матрицу A:\n";
a.Enter_Key();
cout<<"Введите вектор b:\n";
b.Enter_Key();
cout<<"Введите eps:\n";
cin>>eps;
cout<<"Введите начальное приближение x:\n";
xi.Enter_Key();
Matrix m(a + tau*e);

do
{
if (i) xi=xip;
xip = m*xi + tau*b;
i++;

cout<<"-----------------\n";
cout<<"xip:\n";
xip.Print();
cout<<"xi:\n";
xi.Print();
cout<<"abs(xip-xi): "<<abs(xip - xi)<<"\n";
}
while ( abs(xip - xi) > eps );

cout<<"\nОтветы:\n";
xip.Print();
cout<<"\nЧисло итераций: "<<i<<"\n\n";
}

Matrix DOY (Matrix& a)
{
int size = a.GetSize();
Matrix b(size);

b.FillE();
for (int i=0; i<size; i++)
b(i,i) = b(i,i) / a(i,i);

return b;
}

void Yakobi (int size)
{
//xip = (a - d*e)*xi + d*b;
int i = 0;
double eps;

Vector xip(size);
Vector xi(size);
Vector b(size);
Matrix a(size);
Matrix e(size);
Matrix d (size);

cout<<"Метод Якоби\n";
e.FillE();
cout<<"Введите матрицу A:\n";
a.Enter_Key();
cout<<"Введите вектор b:\n";
b.Enter_Key();
cout<<"Введите eps:\n";
cin>>eps;
cout<<"Введите начальное приближение x:\n";
xi.Enter_Key();

d = DOY(a);
Matrix m(a - d*e);

do
{
if (i) xi=xip;
xip = m*xi + d*b;
i++;

cout<<"-----------------\n";
cout<<"xip:\n";
xip.Print();
cout<<"xi:\n";
xi.Print();
cout<<"abs(xip-xi): "<<abs(xip - xi)<<"\n";
}
while ( abs(xip - xi) > eps );

cout<<"\nОтветы:\n";
xip.Print();
cout<<"\nЧисло итераций: "<<i<<"\n\n";
}

int _tmain(int argc, _TCHAR* argv[])
{
int size;
int exit=1;
int method;
char ch;

while (exit)
{
cout<<"Введите число переменных (size): ";
cin>>size;
cout<<"Выберите метод решения:\n1) Метод простых итераций\n2) Метод Якоби\n\n>";
cin>>method;
cout<<"\n\n";

switch (method)
{
case 1:
SIM(size);
break;

case 2:
Yakobi(size);
break;

default:
cout<<"Метод под этим номером не определен!\n";
break;
}

l1: cout<<"Повторить процедуру? [Y/N] ";
cin>>ch;

switch (ch)
{
case 89:
case 121:
break;

case 78:
case 110:
exit = 0;
break;

default:
cout<<"Принимаем ответы только \"Y\" и \"N\"\n";
goto l1;
break;
}
}

return 0;
}

[/codebox]

Файл Classes.h (100% верный, уже не раз использовал его, но для отладки необходим).
[codebox]//Classes.h
#include "stdafx.h"
#include <math.h>

class Vector
{
private:
double *vec;
int size;
public:
Vector(int s);
~Vector();
Vector(const Vector& obj);

Vector& operator=(const Vector& obj);
Vector operator+(const Vector& a);
Vector operator-(const Vector& a);
double& operator[](int i);

void Enter_Key();
void Print();
int GetSize(){return size;}

};

Vector::Vector(int s)
{
if (s>0)
{
vec=new double;
size=s;
}
else
{
if (size<0)
{
vec=new double[abs(s)];
size=abs(s);
}
else
{
vec=new double [1];
size=1;
}
}
}
Vector::~Vector()
{
delete[] vec;
}
Vector::Vector(const Vector& obj)
{
vec=new double[obj.size];
for (int i=0;i<obj.size;i++)
{
vec=obj.vec;
}
size=obj.size;
}
Vector& Vector:perator=(const Vector& obj)
{
delete[] vec;
vec=new double[obj.size];
for (int i=0;i<obj.size;i++)
{
vec=obj.vec;
}
size=obj.size;
return *this;
}
Vector Vector:perator+(const Vector& a)
{
if (size!=a.size) return *this;

Vector b(size);

for (int i=0;i<size;i++)
{
b.vec=vec+a.vec;
}
return b;
}

Vector Vector:perator-(const Vector& a)
{
if (size!=a.size) return *this;

Vector b(size);

for (int i=0;i<size;i++)
{
b.vec=vec-a.vec;
}
return b;
}
double& Vector:perator[](int i)
{
if (i<0) return vec[0];
if (i>=size) return vec[size-1];

return vec;
}
void Vector::Enter_Key()
{
for (int i=0;i<size;i++)
{
cout << "vec[" << i << "]=";
cin >> vec;
}
}
void Vector:rint()
{
for (int i=0;i<size;i++)
{
cout <<vec << " ";
}
cout << "\n";
}

class Matrix
{
private:
double **mat;
int size;
public:
Matrix(int n);
~Matrix();
Matrix(const Matrix& obj);

Matrix& operator=(const Matrix& obj);
Matrix operator+(const Matrix& a);
Matrix operator-(const Matrix& a);
double& operator()(int i,int j);

void Enter_Key();
void Print();
int GetSize(){return size;}
void FillE();
};

Matrix::Matrix(int n)
{
n=abs;

if (n==0) n=1;

mat=new double*[n];
for (int i=0;i<n;i++)
{
mat=new double[n];
}

size=n;
}
Matrix::~Matrix()
{
for (int i=0;i<size;i++)
{
delete[] mat;
}
delete[]mat;
}
Matrix::Matrix(const Matrix& obj)
{
size=obj.size;
mat=new double*[size];
for (int i=0;i<size;i++)
{
mat=new double[size];
}
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
mat[j]=obj.mat[j];
}
}
}
Matrix& Matrix:perator=(const Matrix& obj)
{
for (int i=0;i<size;i++)
{
delete[] mat;
}
delete[]mat;

size=obj.size;
mat=new double*[size];
for (int i=0;i<size;i++)
{
mat=new double[size];
}
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
mat[j]=obj.mat[j];
}
}
return *this;
}
Matrix Matrix:perator+(const Matrix& a)
{
if (size!=a.size) return *this;
Matrix b(size);
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
b.mat[j]=mat[j]+a.mat[j];
}
}
return b;
}
Matrix Matrix:perator-(const Matrix& a)
{
if (size!=a.size) return *this;
Matrix b(size);
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
b.mat[j]=mat[j]-a.mat[j];
}
}
return b;
}
double& Matrix:perator()(int i,int j)
{
if (i<0) i=0;
if (j<0) j=0;
if (i>=size) i=size-1;
if (j>=size) j=size-1;

return *(*(mat+i)+j);
}
void Matrix::Enter_Key()
{
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
cout << "Mat[" << i << "," << j << "]=";
cin >> mat[j];
}
}
}
void Matrix:rint()
{
cout << "Matrix:\n";
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
cout << *(*(mat+i)+j) << " ";
}
cout << "\n";
}
}

void Matrix::FillE()
{
for (int i=0;i<size;i++)
{
for (int j=0;j<size;j++)
{
if (i==j) mat[j]=1;
else mat[j]=0;
}
}
}

Vector operator*(double c,Vector& a)
{
Vector b(a.GetSize());
for (int i=0;i<a.GetSize();i++)
{
b=c*a;
}
return b;
}

Matrix operator*(double c,Matrix& a)
{
Matrix b(a.GetSize());
for (int i=0;i<a.GetSize();i++)
{
for (int j=0;j<a.GetSize();j++)
{
b(i,j)=c*a(i,j);
}
}
return b;
}

Matrix operator*(Matrix& a,Matrix& :blink:
{
if (a.GetSize()!=b.GetSize()) return a;
Matrix c(a.GetSize());
double sum=0;
for (int i=0;i<a.GetSize();i++)
{
for (int j=0;j<a.GetSize();j++)
{
sum=0;
for (int k=0;k<a.GetSize();k++)
{
sum+=a(i,k)*b(k,j);
}
c(i,j)=sum;
}
}
return c;
}
Vector operator*(Matrix& a,Vector&
{
if (a.GetSize()!=b.GetSize()) return b;
Vector c(a.GetSize());
double sum=0;
for (int i=0;i<a.GetSize();i++)
{
sum=0;
for (int j=0;j<a.GetSize();j++)
{
sum+=a(i,j)*b[j];
}
c=sum;
}
return c;
}

double abs(Vector& a)
{
double s=0;
for (int i=0;i<a.GetSize();i++)
{
s+=a*a;
}
return sqrt(s);

}
[/codebox]

chebog · 20.05.2008

сорри, народ. ошибка была в тупой опечатке:

вместо "a - tau*e" надо было писать "e - tau*a".
ну и аналогично для второго метода.

да, еще одно - если кто-то будет это читать...
tau в методе простых итераций в идеале нужно указывать для каждой системы отдельно. обычно его дают в задании. но можно и самим брать. суть такова, что tau должно быть меньше, чем 1/abs(lambda_max), где lambda_max - наибольшее собственное число матрицы А.

Все сервисы Codeby

Поиск

Поиск

Решение СЛАУ методоми простых итераций и Якоби

chebog

chebog