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

Тема в разделе "MS Visual C++", создана пользователем chebog, 20 май 2008.

Статус темы:
Закрыта.
  1. chebog

    chebog Гость

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

    Написал код, который по сути должен решать СЛАУ методами простых итераций и Якоби, но когда пытаюсь решить уравнения, ясно видно, что методы не сходятся... Весь вопрос - почему. По всем условиям должны сходиться. Сейчас отлаживаю, но сдавать уже через 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::eek: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::eek: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::eek: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::eek: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::print()
    {
    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(n);

    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::eek: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[i]=new double[size];
    }
    for (int i=0;i<size;i++)
    {
    for (int j=0;j<size;j++)
    {
    mat[i][j]=obj.mat[i][j];
    }
    }
    return *this;
    }
    Matrix Matrix::eek: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[i][j]=mat[i][j]+a.mat[i][j];
    }
    }
    return b;
    }
    Matrix Matrix::eek: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[i][j]=mat[i][j]-a.mat[i][j];
    }
    }
    return b;
    }
    double& Matrix::eek: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[i][j];
    }
    }
    }
    void Matrix::print()
    {
    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[i][j]=1;
    else mat[i][j]=0;
    }
    }
    }



    Vector operator*(double c,Vector& a)
    {
    Vector b(a.GetSize());
    for (int i=0;i<a.GetSize();i++)
    {
    b[i]=c*a[i];
    }
    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[i]=sum;
    }
    return c;
    }

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

    }
    [/codebox][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i]
     
  2. chebog

    chebog Гость

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

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

    да, еще одно - если кто-то будет это читать...
    tau в методе простых итераций в идеале нужно указывать для каждой системы отдельно. обычно его дают в задании. но можно и самим брать. суть такова, что tau должно быть меньше, чем 1/abs(lambda_max), где lambda_max - наибольшее собственное число матрицы А.
     
Загрузка...
Похожие Темы - Решение СЛАУ методоми
  1. Krex
    Ответов:
    0
    Просмотров:
    1.476
  2. Krex
    Ответов:
    1
    Просмотров:
    1.576
  3. paxac
    Ответов:
    0
    Просмотров:
    45
  4. aameno2
    Ответов:
    0
    Просмотров:
    220
  5. Даниил
    Ответов:
    0
    Просмотров:
    817
Статус темы:
Закрыта.

Поделиться этой страницей