Advertisement
Guest User

Untitled

a guest
Nov 26th, 2014
185
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 15.68 KB | None | 0 0
  1. #define _CRT_SECURE_NO_WARNINGS
  2. #include <iostream>
  3. #include <cstdlib>
  4. #include <cstring>
  5. #include <fstream>
  6. #include <cstdio>
  7. #include <math.h>
  8. using namespace std;
  9.  
  10. class Matrix
  11. {
  12. private:
  13.     int N;
  14.     double **matrix;
  15. public:
  16.     Matrix()
  17.     {
  18.         // реализация конструктора по умолчанию
  19.         N = 0;
  20.         matrix = new double*[N];
  21.         for (int i = 0; i < N; i++)
  22.         {
  23.             matrix[i] = new double[N];
  24.         }
  25.     }
  26.  
  27.     //-----------------------------------------------------------------------------------------------
  28.  
  29.     Matrix(int n)
  30.     {
  31.         // реализация конструктора с параметрами(просто задаем пустую матрицу N*N)
  32.         N = n;
  33.         matrix = new double*[N];
  34.         for (int i = 0; i < N; i++)
  35.         {
  36.             matrix[i] = new double[N];
  37.         }
  38.         for (int i = 0; i < N; i++)
  39.         {
  40.             for (int j = 0; j < N; j++)
  41.             {
  42.                 matrix[i][j] = 0;
  43.             }
  44.         }
  45.     }
  46.  
  47.     //-----------------------------------------------------------------------------------------------
  48.  
  49.     Matrix(const Matrix &nm)
  50.     {
  51.         // реализация конструктора копии
  52.         N = nm.N;
  53.         matrix = new double*[N];
  54.         for (int i = 0; i < N; i++)
  55.         {
  56.             matrix[i] = new double[N];
  57.         }
  58.  
  59.         for (int i = 0; i < N; i++)
  60.         {
  61.             for (int j = 0; j < N; ++j)
  62.             {
  63.                 matrix[i][j] = nm.matrix[i][j];
  64.             }
  65.         }
  66.     }
  67.  
  68.     //-----------------------------------------------------------------------------------------------
  69.  
  70.     ~Matrix()
  71.     {
  72.         // деструктор класса Matrix
  73.         delete[] matrix; // освободить память, удалив матрицу
  74.     }
  75.  
  76.     //-----------------------------------------------------------------------------------------------
  77.  
  78.     Matrix& operator=(const Matrix &nm)
  79.     {
  80.         // присваивание
  81.         N = nm.N;
  82.         for (int i = 0; i < N; i++)
  83.         {
  84.             for (int j = 0; j < N; ++j)
  85.             {
  86.                 matrix[i][j] = nm.matrix[i][j];
  87.             }
  88.         }
  89.         return *this;
  90.     }
  91.  
  92.     //-----------------------------------------------------------------------------------------------
  93.  
  94.     void Matrix::scan()
  95.     {
  96.         //считывание матрицы
  97.         int i = 0, j = 0;
  98.         for (i = 0; i < N; i++)
  99.         {
  100.             for (j = 0; j < N; j++)
  101.             {
  102.                 scanf("%lf", &matrix[i][j]);
  103.             }
  104.         }
  105.     }
  106.  
  107.     //-----------------------------------------------------------------------------------------------
  108.  
  109.     void Matrix::print()
  110.     {
  111.         // вывод матрицы
  112.         int i = 0, j = 0;
  113.         for (i = 0; i < N; i++)
  114.         {
  115.             for (j = 0; j < N; j++)
  116.             {
  117.                 printf("%lf ", matrix[i][j]);
  118.             }
  119.             printf("\n");
  120.         }
  121.     }
  122.  
  123.     //-----------------------------------------------------------------------------------------------
  124.  
  125.     Matrix operator *(const Matrix &nm) const
  126.     {
  127.         //Перемножение матриц(nm слева)
  128.         int i, j, k;
  129.         double dump;
  130.         Matrix retma(N);
  131.         for (i = 0; i < N; i++)
  132.         {
  133.             for (j = 0; j < N; j++)
  134.             {
  135.                 dump = 0;
  136.                 for (k = 0; k < N; k++)
  137.                 {
  138.                     dump += matrix[k][j] * nm.matrix[i][k];
  139.                 }
  140.                 retma.matrix[i][j] = dump;
  141.             }
  142.         }
  143.         return retma;
  144.     }
  145.  
  146.     //-----------------------------------------------------------------------------------------------
  147.  
  148.     Matrix operator *(const double scalar) const
  149.     {
  150.         //произведение матрицы на скаляр
  151.         int i = 0, j = 0;
  152.         double result = 0;
  153.         Matrix retma(N);
  154.         for (i = 0; i < N; i++)
  155.         {
  156.             for (j = 0; j < N; j++)
  157.             {
  158.                 retma.matrix[i][j] = matrix[i][j] * scalar;
  159.             }
  160.  
  161.         }
  162.         return retma;
  163.     }
  164.  
  165.     //-----------------------------------------------------------------------------------------------
  166.  
  167.     Matrix operator +(const Matrix &v2) const
  168.     {
  169.         //Сумма матриц
  170.         Matrix retma(N);
  171.         int i = 0, j = 0;
  172.         for (i = 0; i < N; i++)
  173.         {
  174.             for (j = 0; j < N; j++)
  175.             {
  176.                 retma.matrix[i][j] = matrix[i][j] + v2.matrix[i][j];
  177.             }
  178.         }
  179.         return retma;
  180.     }
  181.  
  182.     //-----------------------------------------------------------------------------------------------
  183.  
  184.     Matrix operator -(const Matrix &v2) const
  185.     {
  186.         //Матрица слева минус матрица справа
  187.         Matrix retma(N);
  188.         int i = 0, j = 0;
  189.         for (i = 0; i < N; i++)
  190.         {
  191.             for (j = 0; j < N; j++)
  192.             {
  193.                 retma.matrix[i][j] = matrix[i][j] - v2.matrix[i][j];
  194.             }
  195.         }
  196.         return retma;
  197.     }
  198.  
  199.     //-----------------------------------------------------------------------------------------------
  200.  
  201.     double Matrix::elem(const int row, int col) const
  202.     {
  203.         //вызов элемента матрицы по индексу(matrix.elem(row, col))
  204.         return matrix[row][col];
  205.     }
  206.  
  207.     //-----------------------------------------------------------------------------------------------
  208.  
  209.     void Matrix::newelem(const int row, int col, double newel)
  210.     {
  211.         //изменение элемента матрицы(matrix.newelem(row, col, new element)
  212.         matrix[row][col] = newel;
  213.     }
  214.  
  215.     //-----------------------------------------------------------------------------------------------
  216.  
  217.     double Matrix::norma1() const
  218.     {
  219.         //Первая матричная норма(максимум из строковых сумм)
  220.         double norma = 0, sum = 0;
  221.         int i = 0, j = 0;
  222.         for (i = 0; i < N; i++)
  223.         {
  224.             for (j = 0; j < N; j++)
  225.             {
  226.                 sum += fabs(matrix[i][j]);
  227.             }
  228.             if (sum>norma)
  229.             {
  230.                 norma = sum;
  231.             }
  232.             sum = 0;
  233.         }
  234.         return norma;
  235.     }
  236.  
  237.     //-----------------------------------------------------------------------------------------------
  238.  
  239.     Matrix Matrix::transp() const
  240.     {
  241.         //Транспонирование матрицы
  242.         int i = 0, j = 0;
  243.         Matrix A = Matrix(N);
  244.         for (i = 0; i < N; i++)
  245.         {
  246.             for (j = 0; j < N; j++)
  247.             {
  248.                 A.matrix[i][j] = matrix[j][i];
  249.             }
  250.         }
  251.         return A;
  252.     }
  253.  
  254.     //-----------------------------------------------------------------------------------------------
  255.  
  256.     double Matrix::norma3() const
  257.     {
  258.         //Первая матричная норма(максимум из столбцевых сумм)
  259.         double norma = 0, sum = 0;
  260.         int i = 0, j = 0;
  261.         for (i = 0; i < N; i++)
  262.         {
  263.             for (j = 0; j < N; j++)
  264.             {
  265.                 sum += fabs(matrix[j][i]);
  266.             }
  267.             if (sum>norma)
  268.             {
  269.                 norma = sum;
  270.             }
  271.             sum = 0;
  272.         }
  273.         return norma;
  274.     }
  275.  
  276.     //-----------------------------------------------------------------------------------------------
  277.  
  278.     Matrix Matrix::rotation(int row, int col)
  279.     {
  280.         //создание матрицы вращения для элемента A(row,col)(matrix.rotation(row, col))
  281.         Matrix retma(N);
  282.         double a = 0, b = 0;
  283.         int i = 0;
  284.         //задание a и b
  285.         a = matrix[col][col] / sqrt(matrix[col][col] * matrix[col][col] + matrix[row][col] * matrix[row][col]);
  286.         b = matrix[row][col] / sqrt(matrix[col][col] * matrix[col][col] + matrix[row][col] * matrix[row][col]);
  287.         //заполнение матрицы единичками
  288.         for (i = 0; i < N; i++)
  289.         {
  290.             retma.matrix[i][i] = 1;
  291.         }
  292.         //запись в матрицу a и b
  293.         retma.matrix[row][row] = a;
  294.         retma.matrix[col][col] = a;
  295.         //условие на элементы(если под диагоналями, то бэшечки распологаются так же, но индексы у нас другие)))))
  296.         retma.matrix[col][row] = b;
  297.         retma.matrix[row][col] = (-1)*b;
  298.  
  299.         return retma;
  300.     }
  301. };
  302.  
  303. //=================================================================================================
  304. //=================================================================================================
  305. //=================================================================================================
  306.  
  307. class Vektor
  308. {
  309. private:
  310.     int N;
  311.     double *vektor;
  312. public:
  313.     Vektor()
  314.     {
  315.         // реализация конструктора по умолчанию
  316.         N = 0;
  317.         vektor = new double[N];
  318.         for (int i = 0; i < N; i++)
  319.         {
  320.             vektor[i] = 0;
  321.         }
  322.     }
  323.  
  324.     //-----------------------------------------------------------------------------------------------
  325.  
  326.     Vektor(int n)
  327.     {
  328.         // реализация конструктора с параметрами(задаем пустой массив размера N)
  329.         N = n;
  330.         vektor = new double[N];
  331.         for (int i = 0; i < N; i++)
  332.         {
  333.             vektor[i] = 0;
  334.         }
  335.     }
  336.  
  337.     //-----------------------------------------------------------------------------------------------
  338.  
  339.     Vektor(const Vektor &nm)
  340.     {
  341.         // реализация конструктора копии
  342.         N = nm.N;
  343.         vektor = new double[N];
  344.         for (int i = 0; i < N; i++)
  345.         {
  346.             vektor[i] = nm.vektor[i];
  347.         }
  348.     }
  349.  
  350.     //-----------------------------------------------------------------------------------------------
  351.  
  352.     ~Vektor()
  353.     {
  354.         // деструктор класса Vektor
  355.         delete[] vektor; // освободить память, удалив вектор
  356.     }
  357.  
  358.     //-----------------------------------------------------------------------------------------------
  359.  
  360.     Vektor& operator=(const Vektor &nm)
  361.     {
  362.         // присваивание
  363.         N = nm.N;
  364.         for (int i = 0; i < N; i++)
  365.         {
  366.             vektor[i] = nm.vektor[i];
  367.         }
  368.         return *this;
  369.     }
  370.  
  371.     //-----------------------------------------------------------------------------------------------
  372.  
  373.     void Vektor::scan()
  374.     {
  375.         //считывание вектора
  376.         for (int i = 0; i < N; i++)
  377.         {
  378.             scanf("%lf", &vektor[i]);
  379.         }
  380.     }
  381.  
  382.     //-----------------------------------------------------------------------------------------------
  383.  
  384.     void Vektor::print()
  385.     {
  386.         // вывод вектора
  387.         for (int i = 0; i < N; i++)
  388.         {
  389.             printf("%lf\n", vektor[i]);
  390.         }
  391.     }
  392.  
  393.     //-----------------------------------------------------------------------------------------------
  394.  
  395.     Vektor operator *(const Matrix &nm) const
  396.     {
  397.         //Перемножение матрицы на вектор(nm слева)
  398.         int i, j;
  399.         double dump;
  400.         Vektor retvek(N);
  401.         for (i = 0; i < N; i++)
  402.         {
  403.             dump = 0;
  404.             for (j = 0; j < N; j++)
  405.             {
  406.                 dump += vektor[j] * nm.elem(i, j);
  407.             }
  408.             retvek.vektor[i] = dump;
  409.         }
  410.         return retvek;
  411.     }
  412.  
  413.     //-----------------------------------------------------------------------------------------------
  414.  
  415.     double operator *(const Vektor &v2) const
  416.     {
  417.         //скалярное произведение векторов
  418.         double result = 0;
  419.         for (int i = 0; i < N; i++)
  420.         {
  421.             result += vektor[i] * v2.vektor[i];
  422.         }
  423.         return result;
  424.     }
  425.  
  426.     //-----------------------------------------------------------------------------------------------
  427.  
  428.     Vektor operator *(const double scalar) const
  429.     {
  430.         //произведение вектора на скаляр
  431.         double result = 0;
  432.         Vektor retvek(N);
  433.         for (int i = 0; i < N; i++)
  434.         {
  435.             retvek.vektor[i] = vektor[i] * scalar;
  436.         }
  437.         return retvek;
  438.     }
  439.  
  440.     //-----------------------------------------------------------------------------------------------
  441.  
  442.     Vektor operator /(const double scalar) const
  443.     {
  444.         //деление вектора на скаляр
  445.         double result = 0;
  446.         Vektor retvek(N);
  447.         for (int i = 0; i < N; i++)
  448.         {
  449.             retvek.vektor[i] = vektor[i] / scalar;
  450.         }
  451.         return retvek;
  452.     }
  453.  
  454.     //-----------------------------------------------------------------------------------------------
  455.  
  456.     Vektor operator +(const Vektor &v2) const
  457.     {
  458.         //Сумма векторов
  459.         Vektor retvek(N);
  460.         for (int i = 0; i < N; i++)
  461.         {
  462.             retvek.vektor[i] = vektor[i] + v2.vektor[i];
  463.         }
  464.         return retvek;
  465.     }
  466.  
  467.     //-----------------------------------------------------------------------------------------------
  468.  
  469.     Vektor operator -(const Vektor &v2) const
  470.     {
  471.         //Вектор слева-вектор справа
  472.         Vektor retvek(N);
  473.         for (int i = 0; i < N; i++)
  474.         {
  475.             retvek.vektor[i] = vektor[i] - v2.vektor[i];
  476.         }
  477.         return retvek;
  478.     }
  479.  
  480.     //-----------------------------------------------------------------------------------------------
  481.  
  482.     double Vektor::norma1() const
  483.     {
  484.         //Первая векторная норма(Сумма модулей координат)
  485.         double norma = 0;
  486.         for (int i = 0; i < N; i++)
  487.         {
  488.             norma += fabs(vektor[i]);
  489.         }
  490.         return norma;
  491.     }
  492.  
  493.     //-----------------------------------------------------------------------------------------------
  494.  
  495.     double Vektor::norma2() const
  496.     {
  497.         //Вторая векторная норма(Корень из суммы квадратов координат)
  498.         double norma = 0;
  499.         for (int i = 0; i < N; i++)
  500.         {
  501.             norma += fabs(vektor[i])*fabs(vektor[i]);
  502.         }
  503.         return sqrt(norma);
  504.     }
  505.  
  506.     //-----------------------------------------------------------------------------------------------
  507.  
  508.     double Vektor::norma3() const
  509.     {
  510.         //Бесконечная норма вектора(максимум из модулей координат)
  511.         double norma = 0;
  512.         int j = 0;
  513.         for (int i = 0; i < N; i++)
  514.         {
  515.             if (norma<fabs(vektor[i]))
  516.             {
  517.                 norma = vektor[i];
  518.                 j = i;
  519.             }
  520.         }
  521.         norma = vektor[j];
  522.         return norma;
  523.     }
  524.  
  525.     //-----------------------------------------------------------------------------------------------
  526.  
  527.     double Vektor::elem(const int row) const
  528.     {
  529.         //вызов элемента векторного столбца по индексу(vektor.elem(row))
  530.         return vektor[row];
  531.     }
  532.  
  533.     //-----------------------------------------------------------------------------------------------
  534.  
  535.     void Vektor::newelem(const int row, double newel)
  536.     {
  537.         //изменение элемента матрицы(vektor.newelem(row, new element)
  538.         vektor[row] = newel;
  539.     }
  540.  
  541.     //-----------------------------------------------------------------------------------------------
  542.  
  543.     Vektor operator ^(const Matrix &nm) const
  544.     {
  545.         //Выражение иксов
  546.         int i, j;
  547.         Vektor retvek(N);
  548.         Matrix newmat = nm;
  549.         //деление на A(i)(j)
  550.         //Столбец свободных коэффициентов
  551.         for (i = 0; i < N; i++)
  552.         {
  553.             vektor[i] = vektor[i] / nm.elem(i, i);
  554.         }
  555.         //матрица
  556.         for (i = 0; i < N; i++)
  557.         {
  558.             for (j = 0; j <= i; j++)
  559.             {
  560.                 newmat.newelem(i, j, newmat.elem(i, j) / newmat.elem(i, i));
  561.             }
  562.         }
  563.         retvek.newelem(0, vektor[0]);//присвоили первый номер
  564.         for (i = 1; i < N; i++)
  565.         {
  566.             retvek.vektor[i] = vektor[i];
  567.             for (j = i - 1; j >= 0; j--)
  568.             {
  569.                 retvek.vektor[i] -= newmat.elem(i, j)*retvek.vektor[j];
  570.             }
  571.         }
  572.         return retvek;
  573.     }
  574. };
  575. double spectr(const Matrix &A, int razm)
  576. {
  577.     double lambdamax = 0;
  578.     int i = 0, q = 0;
  579.     Vektor X = Vektor(razm);
  580.     Vektor XN = Vektor(razm);
  581.     for (i = 0; i < razm; i++)
  582.     {
  583.         X.newelem(i, 1);
  584.     }
  585.     XN = (X*A) / X.norma2();
  586.     while (fabs(XN.norma2() - X.norma2())>0.000001)
  587.     {
  588.         X = XN;
  589.         XN = (X*A) / X.norma2();
  590.         q++;
  591.     }
  592.     lambdamax = XN*XN;
  593.     lambdamax = sqrt(lambdamax);
  594.     return lambdamax;
  595. }
  596. //ПОИСК Q
  597. Matrix findQ(const Matrix &A, int razm)
  598. {
  599.     int e = 0, i = 0, w = 0, q=0;
  600.     Matrix R = A;
  601.     Matrix Q = R.rotation(0, 1);
  602.     Q = Q.transp();
  603.     for (e = 0; e <= 5; e++)//для нормальной точности нужно минимум 6 итераций
  604.     {
  605.         for (q = 0; q < razm; q++)
  606.         {
  607.             for (w = q + 1; w < razm; w++)
  608.             {
  609.                 Matrix SO = R.rotation(q, w); //создание матрицы вращения
  610.                 R = R*SO;//перемножение матриц
  611.                 if ((w != 1) || (q != 0) || (e != 0))
  612.                 {
  613.                     Q = Q*SO.transp();
  614.                 }
  615.             }
  616.         }
  617.     }
  618.     return Q;
  619. }
  620. int main()
  621. {
  622.     freopen("in.txt", "r", stdin);
  623.     freopen("out.txt", "w", stdout);
  624.     int razm = 0, i=0, j=0, k=0, h=0;
  625.     double maxnorm = 0;
  626.     scanf("%d", &razm);
  627.     Matrix A = Matrix(razm);
  628.     A.scan();
  629.     Vektor B = Vektor(razm);
  630.     Matrix Q = Matrix(razm);
  631.     do
  632.     {
  633.         maxnorm = 0;
  634.         for (j = 0; j < razm-1; j++)
  635.         {
  636.             for (i = razm-1; i > j; i--)
  637.             {
  638.                 B.newelem(i, A.elem(i, j));
  639.             }
  640.             if (B.norma2()>maxnorm)
  641.             {
  642.                 maxnorm = B.norma2();
  643.             }
  644.             for (k = 0; k < razm - 1; k++)
  645.             {
  646.                 B.newelem(k, 0);
  647.             }
  648.         }
  649.         Q = findQ(A, razm);
  650.         A=Q*(A*Q.transp());
  651.         h++;
  652.     } while (maxnorm>0.1);
  653.     A.print();
  654.     printf("\n");
  655.     printf("%d", h);
  656.  
  657.     return 0;
  658. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement