Фрагмент для ознакомления
n}for (int i = 0; i < n; i++){for (int j = 0; j < n; j++){B[i][j] = A[i][j];}for (int j = n; j < 2 * n; j++){if (i == j - n) B[i][j] = 1;}}for (int i = n - 1; i > 0; i--){if (B[i - 1, 0]>= B[i, 0]) continue;for (int j = 0; j < 2 * n; j++){float temp = B[i][j];B[i][j] = B[i - 1][j];B[i - 1][j] = temp;}}for (int i = 0; i < n; i++){for (int j = 0; j < n; j++){if (i == j) continue;float temp = B[j][i] / B[i][i];for (int k = 0; k < 2 * n; k++){B[j][k] -= B[i][k] * temp;}}}for (int i = 0; i < n; i++){float temp = B[i][i];for (int j = 0; j < 2 * n; j++){B[i][j] = B[i][j] / temp;}}auto C = std::make_shared< std::shared_ptr[] >(n);for (size_t i = 0; i < n; i++) {C[i]= std::make_shared(n);}for (int i = 0; i < n; i++){for (int j = n; j < 2 * n; j++){C[i][j - n] = B[i][j];}}return C;}public: static std::shared_ptr[]> Sqrt(intn, std::shared_ptr[]> A){std::shared_ptr[]> B = Duplicate(n, A);for (int i = 0; i < n; i++){if (A[i][i] <= 0) { thrownew std::exception("Divide by zero"); };B[i][i] = (float)std::sqrt(A[i][i]);}return B;}public: static std::string toString(introwCount, intcolumnCount, std::shared_ptr[]> A){std::string text = "";for (int i = 0; i < rowCount; i++){if (i > 0) text +=',';text +='{';for (int j = 0; j < columnCount; j++){if (j > 0) text +=',';text +=A[i][j];}text +='}';}return text;}// распечаткавектораpublic: static std::string toString(intn, std::shared_ptr a){std::string text = "{";for (int i = 0; i < n; i++){if (i > 0) text +=',';text +=a[i];}text +='}';return text;}};SolveEigenvaluesEigenvectors.h#pragmaonce#include#include#include#include"LinearAlgebra.h"classSolveEigenvaluesEigenvectors{// A - исходнаяматрица// iterations - числоитераций// eigenvalues - собственные значения// eigenvectors - собственные вектораpublic: staticvoid find(intn,std::shared_ptr[]> A,intiterations,std::shared_ptr[]>& eigenvalues,std::shared_ptr[]>& eigenvectors){std::shared_ptr[]> B = LinearAlgebra::Duplicate(n, A);std::shared_ptr[]> U = LinearAlgebra::Identity(n);for (int i = 0; i < iterations; i++){std::shared_ptr[]> Q, R;decomposition(n, B, Q, R);B =LinearAlgebra::Product(n, R, Q);U =LinearAlgebra::Product(n, U, Q);}eigenvalues= std::make_shared< std::shared_ptr[] >(n);for (size_t i = 0; i < n; i++) {eigenvalues[i]= std::make_shared(n);}for (int i = 0; i < n; i++){eigenvalues[i][i] = B[i][i];}eigenvectors= U;}private: staticvoid decomposition(intn,std::shared_ptr[]> A,std::shared_ptr[]>& Q,std::shared_ptr[]>& R){std::shared_ptr[]> U = LinearAlgebra::Duplicate(n, A);for (int j = 1; j < n; j++){std::shared_ptr u = LinearAlgebra::GetColumn(n, U, j);std::shared_ptr v = LinearAlgebra::GetColumn(n, U, j);for (int k = j - 1; k >= 0; k--){std::shared_ptr uk = LinearAlgebra::GetColumn(n, U, k);u =LinearAlgebra::Subtract(n, u, LinearAlgebra::Project(n, uk, v));}for (int i = 0; i < n; i++){U[i][j] = u[i];}}for (int j = 0; j < n; j++){std::shared_ptr u = LinearAlgebra::GetColumn(n, U, j);float magnitude = LinearAlgebra::Magnitude(n, u);for (int i = 0; i < n; i++){U[i][j] = u[i] / magnitude;}}Q= U;R=LinearAlgebra::Product(n, LinearAlgebra::Transpose(n, Q), A);}};JacobSolve.h#pragmaonce#include#include"Matrix.h"#include"LinearAlgebra.h"classJacobSolve{private:staticbool symmetrycheck(Matrix& m) {for (int i = 0; i < m.n - 1; i++) {for (int j = i + 1; j < m.n; j++) {if (m.mat[i][j] != m.mat[j][i]) {returnfalse;}}}returntrue;}public:static std::string solve(Matrix& m) {if (!symmetrycheck(m)) {std::cout <<"матрицанесимметрична"<< std::endl;return 0;}auto Delta = std::make_shared(m.n);for (int minor_i = 1; minor_i <= m.n; minor_i++) {auto submat = std::make_shared< std::shared_ptr[] >(minor_i);for (size_t i = 0; i < minor_i; i++) {submat[i]= std::make_shared(minor_i);}for (int i = 0; i < minor_i; i++) {for (int j = 0; j < minor_i; j++) {submat[i][j] = m.mat[i][j];}}float det = LinearAlgebra::Determinant(minor_i, submat);Delta[minor_i - 1] = det;}std::string s = "Q(";for (int i = 0; i < m.n; i++) {s +="y"+ std::to_string(i + 1);if (i != m.n - 1) s +=",";}s +=")=";for (int i = 0; i < m.n; i++) {float koeff = 0;if (i > 0) {koeff = Delta[i] / Delta[i - 1];}else {koeff = Delta[0];}if (koeff >= 0)s +=" + ";elses +=" - ";s += std::to_string(abs(koeff)) +"*y"+ std::to_string(i + 1) +"^2";}return s;}};main.cpp#include#include"LinearAlgebra.h"#include"SolveEigenvaluesEigenvectors.h"#include"Matrix.h"#include"JacobSolve.h"#include#includeint main(){std::system("chcp 1251");int n = 3;auto matrix = std::make_shared< std::shared_ptr[] >(n);for (size_t i = 0; i < n; i++) {matrix[i]= std::make_shared(n);}Matrix mat1(n, matrix);std::cout <<"Введитематрицувформате:"<< std::endl <<"[размерность]"<< std::endl <<"[a11] ... [a1n]"<< std::endl<<"..."<< std::endl <<"[an1] ... [ann]"<< std::endl;std::cin >> mat1;//float matrix_data[3][3] = {//{1, 2, 3},//{2, 3, 5},//{3, 5, 6}//};//for (size_t i = 0; i < n; i++) {//for (size_t j = 0; j < n; j++) {//matrix[i][j] = matrix_data[i][j];//}//}std::cout <<"Нахождение квадратичной формы, метод Якоба:"<< std::endl;JacobSolve jacobSolve;std::string jbStr = jacobSolve.solve(mat1);std::cout <<"Каноническая запись:"<< std::endl;std::cout << jbStr;std::cout << std::endl;std::cout << std::endl;// число итерацийint numberOfIterations = 10;std::shared_ptr[]> eigenvalues;std::shared_ptr[]> eigenvectors;// выполняем поиск собственных веторов и чиселSolveEigenvaluesEigenvectors::find(mat1.n,mat1.mat,numberOfIterations,eigenvalues,eigenvectors);std::cout <<"Тензор"<< std::endl;std::cout << mat1 << std::endl;for (size_t i = 0; i < n; i++) {for (size_t j = 0; j < n; j++) {std::cout << matrix[i][j]<<" ";}std::cout << std::endl;}std::cout << std::endl;std::cout << std::endl;auto sortEV = std::vector< std::pair >();for (int i = 0; i < n; i++) {sortEV.push_back(std::pair(i, eigenvalues[i][i]));}std::sort(sortEV.begin(), sortEV.end(), [](const std::pair& x,const std::pair& y){return (x.second > y.second);});std::cout <<"Собственные значения (отсортированы по убыванию)"<< std::endl;/*for (size_t i = 0; i < n; i++) {for (size_t j = 0; j < n; j++) {std::cout << eigenvalues[i][j] << " ";}std::cout << std::endl;}*/for (auto it = sortEV.begin(); it != sortEV.end(); it++) {int i = it->first;std::cout << eigenvalues[i][i]<<" ";}std::cout << std::endl;std::cout << std::endl;std::cout <<"Собственные вектора"<< std::endl;//for (size_t i = 0; i < n; i++) {//for (size_t j = 0; j < n; j++) {//std::cout << eigenvectors[i][j] << " ";//}//std::cout << std::endl;//}for (auto it = sortEV.begin(); it != sortEV.end(); it++) {int i = it->first;for (size_t j = 0; j < n; j++) {std::cout << eigenvectors[j][i]<<" ";}std::cout << std::endl;}std::cout << std::endl;std::cout <<"Нахождение квадратичной формы, собственные значения и вектора:"<< std::endl;for (size_t j = 0; j < n; j++) {float val = 0;for (size_t i = 0; i < n; i++) {val += eigenvectors[i][j] * eigenvectors[i][j];}val = sqrt(val);for (size_t i = 0; i < n; i++) {eigenvectors[i][j] /= val;}}std::cout <<"Каноническаязапись:"<< std::endl;std::cout <<"Q(";for (int i = 0; i < n; i++) {std::cout <<"y"<< (i + 1);if (i != n - 1) std::cout <<",";}std::cout <<")=";for (int i = 0; i < n; i++) {if (eigenvalues[i][i] >= 0)std::cout <<" + ";elsestd::cout <<" - ";std::cout << abs(eigenvalues[i][i]) <<"*y"<< (i + 1) <<"^2";}std::cout << std::endl;std::cout <<"P=";Matrix mNormEV(n, eigenvectors);std::cout << mNormEV;std::cout << std::endl;}
