Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #define _CRT_SECURE_NO_WARNINGS
- #include <iostream>
- #include <cstdlib>
- #include <cstring>
- #include <fstream>
- #include <cstdio>
- #include <math.h>
- using namespace std;
- class Matrix
- {
- private:
- int N;
- double **matrix;
- public:
- Matrix()
- {
- // реализация конструктора по умолчанию
- N = 0;
- matrix = new double*[N];
- for (int i = 0; i < N; i++)
- {
- matrix[i] = new double[N];
- }
- }
- //-----------------------------------------------------------------------------------------------
- Matrix(int n)
- {
- // реализация конструктора с параметрами(просто задаем пустую матрицу N*N)
- N = n;
- matrix = new double*[N];
- for (int i = 0; i < N; i++)
- {
- matrix[i] = new double[N];
- }
- for (int i = 0; i < N; i++)
- {
- for (int j = 0; j < N; j++)
- {
- matrix[i][j] = 0;
- }
- }
- }
- //-----------------------------------------------------------------------------------------------
- Matrix(const Matrix &nm)
- {
- // реализация конструктора копии
- N = nm.N;
- matrix = new double*[N];
- for (int i = 0; i < N; i++)
- {
- matrix[i] = new double[N];
- }
- for (int i = 0; i < N; i++)
- {
- for (int j = 0; j < N; ++j)
- {
- matrix[i][j] = nm.matrix[i][j];
- }
- }
- }
- //-----------------------------------------------------------------------------------------------
- ~Matrix()
- {
- // деструктор класса Matrix
- delete[] matrix; // освободить память, удалив матрицу
- }
- //-----------------------------------------------------------------------------------------------
- Matrix& operator=(const Matrix &nm)
- {
- // присваивание
- N = nm.N;
- for (int i = 0; i < N; i++)
- {
- for (int j = 0; j < N; ++j)
- {
- matrix[i][j] = nm.matrix[i][j];
- }
- }
- return *this;
- }
- //-----------------------------------------------------------------------------------------------
- void Matrix::scan()
- {
- //считывание матрицы
- int i = 0, j = 0;
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- scanf("%lf", &matrix[i][j]);
- }
- }
- }
- //-----------------------------------------------------------------------------------------------
- void Matrix::print()
- {
- // вывод матрицы
- int i = 0, j = 0;
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- printf("%lf ", matrix[i][j]);
- }
- printf("\n");
- }
- }
- //-----------------------------------------------------------------------------------------------
- Matrix operator *(const Matrix &nm) const
- {
- //Перемножение матриц(nm слева)
- int i, j, k;
- double dump;
- Matrix retma(N);
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- dump = 0;
- for (k = 0; k < N; k++)
- {
- dump += matrix[k][j] * nm.matrix[i][k];
- }
- retma.matrix[i][j] = dump;
- }
- }
- return retma;
- }
- //-----------------------------------------------------------------------------------------------
- Matrix operator *(const double scalar) const
- {
- //произведение матрицы на скаляр
- int i = 0, j = 0;
- double result = 0;
- Matrix retma(N);
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- retma.matrix[i][j] = matrix[i][j] * scalar;
- }
- }
- return retma;
- }
- //-----------------------------------------------------------------------------------------------
- Matrix operator +(const Matrix &v2) const
- {
- //Сумма матриц
- Matrix retma(N);
- int i = 0, j = 0;
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- retma.matrix[i][j] = matrix[i][j] + v2.matrix[i][j];
- }
- }
- return retma;
- }
- //-----------------------------------------------------------------------------------------------
- Matrix operator -(const Matrix &v2) const
- {
- //Матрица слева минус матрица справа
- Matrix retma(N);
- int i = 0, j = 0;
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- retma.matrix[i][j] = matrix[i][j] - v2.matrix[i][j];
- }
- }
- return retma;
- }
- //-----------------------------------------------------------------------------------------------
- double Matrix::elem(const int row, int col) const
- {
- //вызов элемента матрицы по индексу(matrix.elem(row, col))
- return matrix[row][col];
- }
- //-----------------------------------------------------------------------------------------------
- void Matrix::newelem(const int row, int col, double newel)
- {
- //изменение элемента матрицы(matrix.newelem(row, col, new element)
- matrix[row][col] = newel;
- }
- //-----------------------------------------------------------------------------------------------
- double Matrix::norma1() const
- {
- //Первая матричная норма(максимум из строковых сумм)
- double norma = 0, sum = 0;
- int i = 0, j = 0;
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- sum += fabs(matrix[i][j]);
- }
- if (sum>norma)
- {
- norma = sum;
- }
- sum = 0;
- }
- return norma;
- }
- //-----------------------------------------------------------------------------------------------
- Matrix Matrix::transp() const
- {
- //Транспонирование матрицы
- int i = 0, j = 0;
- Matrix A = Matrix(N);
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- A.matrix[i][j] = matrix[j][i];
- }
- }
- return A;
- }
- //-----------------------------------------------------------------------------------------------
- double Matrix::norma3() const
- {
- //Первая матричная норма(максимум из столбцевых сумм)
- double norma = 0, sum = 0;
- int i = 0, j = 0;
- for (i = 0; i < N; i++)
- {
- for (j = 0; j < N; j++)
- {
- sum += fabs(matrix[j][i]);
- }
- if (sum>norma)
- {
- norma = sum;
- }
- sum = 0;
- }
- return norma;
- }
- //-----------------------------------------------------------------------------------------------
- Matrix Matrix::rotation(int row, int col)
- {
- //создание матрицы вращения для элемента A(row,col)(matrix.rotation(row, col))
- Matrix retma(N);
- double a = 0, b = 0;
- int i = 0;
- //задание a и b
- a = matrix[col][col] / sqrt(matrix[col][col] * matrix[col][col] + matrix[row][col] * matrix[row][col]);
- b = matrix[row][col] / sqrt(matrix[col][col] * matrix[col][col] + matrix[row][col] * matrix[row][col]);
- //заполнение матрицы единичками
- for (i = 0; i < N; i++)
- {
- retma.matrix[i][i] = 1;
- }
- //запись в матрицу a и b
- retma.matrix[row][row] = a;
- retma.matrix[col][col] = a;
- //условие на элементы(если под диагоналями, то бэшечки распологаются так же, но индексы у нас другие)))))
- retma.matrix[col][row] = b;
- retma.matrix[row][col] = (-1)*b;
- return retma;
- }
- };
- //=================================================================================================
- //=================================================================================================
- //=================================================================================================
- class Vektor
- {
- private:
- int N;
- double *vektor;
- public:
- Vektor()
- {
- // реализация конструктора по умолчанию
- N = 0;
- vektor = new double[N];
- for (int i = 0; i < N; i++)
- {
- vektor[i] = 0;
- }
- }
- //-----------------------------------------------------------------------------------------------
- Vektor(int n)
- {
- // реализация конструктора с параметрами(задаем пустой массив размера N)
- N = n;
- vektor = new double[N];
- for (int i = 0; i < N; i++)
- {
- vektor[i] = 0;
- }
- }
- //-----------------------------------------------------------------------------------------------
- Vektor(const Vektor &nm)
- {
- // реализация конструктора копии
- N = nm.N;
- vektor = new double[N];
- for (int i = 0; i < N; i++)
- {
- vektor[i] = nm.vektor[i];
- }
- }
- //-----------------------------------------------------------------------------------------------
- ~Vektor()
- {
- // деструктор класса Vektor
- delete[] vektor; // освободить память, удалив вектор
- }
- //-----------------------------------------------------------------------------------------------
- Vektor& operator=(const Vektor &nm)
- {
- // присваивание
- N = nm.N;
- for (int i = 0; i < N; i++)
- {
- vektor[i] = nm.vektor[i];
- }
- return *this;
- }
- //-----------------------------------------------------------------------------------------------
- void Vektor::scan()
- {
- //считывание вектора
- for (int i = 0; i < N; i++)
- {
- scanf("%lf", &vektor[i]);
- }
- }
- //-----------------------------------------------------------------------------------------------
- void Vektor::print()
- {
- // вывод вектора
- for (int i = 0; i < N; i++)
- {
- printf("%lf\n", vektor[i]);
- }
- }
- //-----------------------------------------------------------------------------------------------
- Vektor operator *(const Matrix &nm) const
- {
- //Перемножение матрицы на вектор(nm слева)
- int i, j;
- double dump;
- Vektor retvek(N);
- for (i = 0; i < N; i++)
- {
- dump = 0;
- for (j = 0; j < N; j++)
- {
- dump += vektor[j] * nm.elem(i, j);
- }
- retvek.vektor[i] = dump;
- }
- return retvek;
- }
- //-----------------------------------------------------------------------------------------------
- double operator *(const Vektor &v2) const
- {
- //скалярное произведение векторов
- double result = 0;
- for (int i = 0; i < N; i++)
- {
- result += vektor[i] * v2.vektor[i];
- }
- return result;
- }
- //-----------------------------------------------------------------------------------------------
- Vektor operator *(const double scalar) const
- {
- //произведение вектора на скаляр
- double result = 0;
- Vektor retvek(N);
- for (int i = 0; i < N; i++)
- {
- retvek.vektor[i] = vektor[i] * scalar;
- }
- return retvek;
- }
- //-----------------------------------------------------------------------------------------------
- Vektor operator /(const double scalar) const
- {
- //деление вектора на скаляр
- double result = 0;
- Vektor retvek(N);
- for (int i = 0; i < N; i++)
- {
- retvek.vektor[i] = vektor[i] / scalar;
- }
- return retvek;
- }
- //-----------------------------------------------------------------------------------------------
- Vektor operator +(const Vektor &v2) const
- {
- //Сумма векторов
- Vektor retvek(N);
- for (int i = 0; i < N; i++)
- {
- retvek.vektor[i] = vektor[i] + v2.vektor[i];
- }
- return retvek;
- }
- //-----------------------------------------------------------------------------------------------
- Vektor operator -(const Vektor &v2) const
- {
- //Вектор слева-вектор справа
- Vektor retvek(N);
- for (int i = 0; i < N; i++)
- {
- retvek.vektor[i] = vektor[i] - v2.vektor[i];
- }
- return retvek;
- }
- //-----------------------------------------------------------------------------------------------
- double Vektor::norma1() const
- {
- //Первая векторная норма(Сумма модулей координат)
- double norma = 0;
- for (int i = 0; i < N; i++)
- {
- norma += fabs(vektor[i]);
- }
- return norma;
- }
- //-----------------------------------------------------------------------------------------------
- double Vektor::norma2() const
- {
- //Вторая векторная норма(Корень из суммы квадратов координат)
- double norma = 0;
- for (int i = 0; i < N; i++)
- {
- norma += fabs(vektor[i])*fabs(vektor[i]);
- }
- return sqrt(norma);
- }
- //-----------------------------------------------------------------------------------------------
- double Vektor::norma3() const
- {
- //Бесконечная норма вектора(максимум из модулей координат)
- double norma = 0;
- int j = 0;
- for (int i = 0; i < N; i++)
- {
- if (norma<fabs(vektor[i]))
- {
- norma = vektor[i];
- j = i;
- }
- }
- norma = vektor[j];
- return norma;
- }
- //-----------------------------------------------------------------------------------------------
- double Vektor::elem(const int row) const
- {
- //вызов элемента векторного столбца по индексу(vektor.elem(row))
- return vektor[row];
- }
- //-----------------------------------------------------------------------------------------------
- void Vektor::newelem(const int row, double newel)
- {
- //изменение элемента матрицы(vektor.newelem(row, new element)
- vektor[row] = newel;
- }
- //-----------------------------------------------------------------------------------------------
- Vektor operator ^(const Matrix &nm) const
- {
- //Выражение иксов
- int i, j;
- Vektor retvek(N);
- Matrix newmat = nm;
- //деление на A(i)(j)
- //Столбец свободных коэффициентов
- for (i = 0; i < N; i++)
- {
- vektor[i] = vektor[i] / nm.elem(i, i);
- }
- //матрица
- for (i = 0; i < N; i++)
- {
- for (j = 0; j <= i; j++)
- {
- newmat.newelem(i, j, newmat.elem(i, j) / newmat.elem(i, i));
- }
- }
- retvek.newelem(0, vektor[0]);//присвоили первый номер
- for (i = 1; i < N; i++)
- {
- retvek.vektor[i] = vektor[i];
- for (j = i - 1; j >= 0; j--)
- {
- retvek.vektor[i] -= newmat.elem(i, j)*retvek.vektor[j];
- }
- }
- return retvek;
- }
- };
- double spectr(const Matrix &A, int razm)
- {
- double lambdamax = 0;
- int i = 0, q = 0;
- Vektor X = Vektor(razm);
- Vektor XN = Vektor(razm);
- for (i = 0; i < razm; i++)
- {
- X.newelem(i, 1);
- }
- XN = (X*A) / X.norma2();
- while (fabs(XN.norma2() - X.norma2())>0.000001)
- {
- X = XN;
- XN = (X*A) / X.norma2();
- q++;
- }
- lambdamax = XN*XN;
- lambdamax = sqrt(lambdamax);
- return lambdamax;
- }
- //ПОИСК Q
- Matrix findQ(const Matrix &A, int razm)
- {
- int e = 0, i = 0, w = 0, q=0;
- Matrix R = A;
- Matrix Q = R.rotation(0, 1);
- Q = Q.transp();
- for (e = 0; e <= 5; e++)//для нормальной точности нужно минимум 6 итераций
- {
- for (q = 0; q < razm; q++)
- {
- for (w = q + 1; w < razm; w++)
- {
- Matrix SO = R.rotation(q, w); //создание матрицы вращения
- R = R*SO;//перемножение матриц
- if ((w != 1) || (q != 0) || (e != 0))
- {
- Q = Q*SO.transp();
- }
- }
- }
- }
- return Q;
- }
- int main()
- {
- freopen("in.txt", "r", stdin);
- freopen("out.txt", "w", stdout);
- int razm = 0, i=0, j=0, k=0, h=0;
- double maxnorm = 0;
- scanf("%d", &razm);
- Matrix A = Matrix(razm);
- A.scan();
- Vektor B = Vektor(razm);
- Matrix Q = Matrix(razm);
- do
- {
- maxnorm = 0;
- for (j = 0; j < razm-1; j++)
- {
- for (i = razm-1; i > j; i--)
- {
- B.newelem(i, A.elem(i, j));
- }
- if (B.norma2()>maxnorm)
- {
- maxnorm = B.norma2();
- }
- for (k = 0; k < razm - 1; k++)
- {
- B.newelem(k, 0);
- }
- }
- Q = findQ(A, razm);
- A=Q*(A*Q.transp());
- h++;
- } while (maxnorm>0.1);
- A.print();
- printf("\n");
- printf("%d", h);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement