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

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

Форум программистов > Центр помощи > [C++]Решение СЛАУ методом простых итераций (Якоби)

Автор: Alleut 13.5.2008, 02:11

http://ru.wikipedia.org/wiki/Метод_Якоби - описание метода

Задача:
Решить систему линейных уравнений Ax = b методом простых итераций, приведя ее к виду , x = Cx + d,
если A = D + kG, k = 0...5

Код


/ 1.0  -0.2   0.1\
|-0.1   1.0  -0.1| = D
\-0.3   0.2  -1.0/

/0.1  0      0   \
|0     0.1   0    | = G
\0     0     0.1 /

/0.4\
|0.8|  = b
\0.2/

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

Добавлено через 3 минуты и 31 секунду
Вот набросок:

Код


class SqMatrix {
protected:
    double** data;
    int size;
public:
    SqMatrix();
    SqMatrix(int); //constructor
    double element(int, int);
    void setElement(int, int, double);    
    ~SqMatrix(); //destructor    
    SqMatrix operator+(SqMatrix);
};

SqMatrix::SqMatrix(){
    size = 1;
    data = new double*[size];
    for (int i = 0; i<size; i++){
      data[i] = new double[size];
      for(int j = 0; j<size; j++)
        data[i][j] = 0.0;
    };
};

SqMatrix::SqMatrix(int s){   //constructor real
    size = s;
    data = new double*[s];
    for (int i = 0; i<s; i++){
      data[i] = new double[s];
      for(int j = 0; j < s; j++)
        data[i][j] = 0.0;
    };
};

double SqMatrix::element(int i, int j){  //??
  return data[i][j];
};

void SqMatrix::setElement(int i, int j, double value){
    data[i][j] = value;
};

SqMatrix::~SqMatrix(){  //деструктор глючит
    for(int i = 0; i < size; i++)
        delete[] data[i];   //глючит где-то здесь
    delete[] data;                  //глючит где-то здесь
    };

SqMatrix SqMatrix::operator+ (SqMatrix B){     //насчет этой хрени я тоже не уверен :)
    SqMatrix C;
    for (int i = 0; i<B.size; i++){
        for (int j = 0; j<B.size; j++){
            C.data[i][j] = data[i][j] + B.data[i][j];
        };
    };
    return C;
};

Автор: Alleut 13.5.2008, 22:19

Упростил класс, повысив его кривизну. Добавил еще класс векторов.
Нужно теперь перегрузить умножения для матриц, и собственно реализовать сам алгоритм через это умножение, но я не очень понимаю критерий остановки рекурсии...

Код


class SqMat {
public:
    SqMat();
    double data[3][3];
    double element(int, int);
    void setElement(int, int, double);    
    SqMat NumMultMat(int, SqMat);
    SqMat operator+(SqMat);
};

SqMat::SqMat(){
    for(int i = 0; i<3; i++){
        for(int j = 0; j<3; j++){
        data[i][j] = 0.0;
        };
    };
};

double SqMat::element(int i, int j){
  return data[i][j];
};

void SqMat::setElement(int i, int j, double value){
    data[i][j] = value;
};

SqMat SqMat::NumMultMat(int k, SqMat A){
    std::cout << k;
    SqMat temp;
    for (int i = 0; i < 3; i++){
        for (int j = 0; j < 3; j++){
            temp.data[i][j] = k*A.data[i][j];
        };
    };
    for (int i = 0; i < 3; i++){
        for (int j = 0; j < 3; j++){
            A.data[i][j] = temp.data[i][j];
        };
    };
    return A;
};

SqMat SqMat::operator+ (SqMat B){
    SqMat C;
    for (int i = 0; i<3; i++){
        for (int j = 0; j<3; j++){
            C.data[i][j] = data[i][j] + B.data[i][j];
        };
    };
    return C;
};

class Vect {
public:
    Vect ();
    double Vdata[3][1];
    double Velement(int);
    void VsetElement(int, double);
};

Vect::Vect(){
    for(int i = 0; i<3; i++){
        Vdata[i][1] = 0.0;
    };
};

double Vect::Velement(int i){
  return Vdata[i][1];
};

void Vect::VsetElement(int i, double value){
    Vdata[i][1] = value;
};

Автор: Alleut 17.5.2008, 12:54

подскажите плиз какой будет критерий остановки рекурсии x = Cx + d

Автор: Alleut 21.5.2008, 22:22

вот собственно что получится, если писать топорно:
юзайте кому надо, единственное, фиговина работает для заданной в коде матрицы, поэтому можете
1) менять матрицу в коде
2) написать цикл для ввода значений
ПС: нет механизма для несходящихся решений. сами придумывайте если понадобится

Код


#include "stdafx.h"
#include "iostream"
#include "math.h"
#include "conio.h"
using namespace std;
const int n = 3;
const double eps = 0.000000000000001;
 
void main(){
    double A[n][n];        //A = D + k*G
    double D[n][n] = {{1.0,-0.2,0.1},{-0.1,1.0,-0.1},{-0.3,0.2,-1.0}}; 
    double G[n][n] = {{0.1,0.0,0.0},{0.0,0.1,0.0},{0.0,0.0,0.1}};
    double C[n][n];
    double d[n];
    double b[n] = {0.4,0.8,0.2};
    double xk[n];
    double xk1[n];
    double norm[n];
    double tmp;
    double norm_c = 0;
    double max = 0;
    int k;
    double test[n];
    
    cout << "Enter k:";
    cin >> k;

    for(int i = 0; i < n; i++){
        for(int j = 0; j < n;j++){
            A[i][j] = D[i][j]+k*G[i][j];
        };
    };
 
    for(int k = 0; k < n; k++){
        d[k] = b[k]/A[k][k];
    };

    for(int i = 0; i < n; i++){
        for(int j = 0; j < n;j++){
            C[i][j] = (-1)*A[i][j]/A[i][i]; 
        };
        C[i][i] = 0;
    };

    for(int i =  0; i < n;i++){
        
        for(int j = 0; j < n; j++){
            norm_c = fabs(C[i][j]);
        };
    };

    for(int i =  0; i < n;i++){
        xk[i] = d[i];
    };

    while(true){
        for(int i =  0; i < n;i++){
            xk1[i] = d[i];
            tmp = 0;
            for(int j = 0; j < n; j++){
                tmp += xk[j] * C[i][j];
            };
        xk1[i] += tmp;
        };

        for(int i = 0; i < n; i++){
            norm[i] = xk1[i] - xk[i];
        };
        tmp = 0;

        for(int i = 0; i < n; i++){
            if (tmp < fabs(norm[i]))
                tmp = fabs(norm[i]);
        };
        
        if(norm_c < 1)
            break;

        for(int i = 0; i < n; i++){
            xk[i] = xk1[i];
        };
    };

    for(int i = 0; i < n; i++){
        cout<< "x[" << i+1<< "]: " << xk1[i]<< endl;
    };
    
        
_getch();
};

Добавлено через 33 секунды
Мог попасться мусор в коде