API tools faq
paste
Advertisement
egraPA

Untitled

egraPA
Apr 30th, 2024
682
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 17.99 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <vector>
  3. #include <math.h>
  4. #include <complex>
  5. #include <stdlib.h>
  6. #include <sstream>
  7. #include <iomanip>
  8.  
  9. using namespace std;
  10. class Vector;
  11.  
  12. class Matrix {
  13.     protected:
  14.         vector<vector<double>> data;
  15.         int n;
  16.         int m;
  17.     public:
  18.         Matrix() {
  19.             n = 0;
  20.             m = 0;
  21.         }
  22.         Matrix (int x, int y) {
  23.             n = x;
  24.             m = y;
  25.             data.assign(n, vector<double>(m, 0));
  26.         }
  27.         Matrix& operator+ (const Matrix& matrix) {
  28.             if (matrix.m != m || matrix.n != n) {
  29.                 throw runtime_error("Error: the dimensional problem occurred");
  30.             }
  31.             Matrix *ans = new Matrix(n, m);
  32.             for (int i = 0; i < n; i++) {
  33.                 for (int j = 0; j < m; j++) {
  34.                     ans->data[i][j] = data[i][j] + matrix.data[i][j];
  35.                 }
  36.             }
  37.             return *ans;
  38.         }
  39.         Matrix& operator- (const Matrix& matrix) {
  40.             if (matrix.m != m || matrix.n != n) {
  41.                 throw runtime_error("Error: the dimensional problem occurred");
  42.             }
  43.             Matrix *ans = new Matrix(n, m);
  44.             for (int i = 0; i < n; i++) {
  45.                 for (int j = 0; j < m; j++) {
  46.                     ans->data[i][j] = data[i][j] - matrix.data[i][j];
  47.                 }
  48.             }
  49.             return *ans;
  50.         }
  51.         Matrix& operator* (const Matrix& matrix) {
  52.             if (matrix.n != m) {
  53.                 throw runtime_error("Error: the dimensional problem occurred");
  54.             }
  55.             Matrix *ans = new Matrix(n, matrix.m);
  56.             for (int i = 0; i < n; i++) {
  57.                 for (int j = 0; j < matrix.m; j++) {
  58.                     for (int k = 0; k < m; k++) {
  59.                         ans->data[i][j] += data[i][k]*matrix.data[k][j];
  60.                     }
  61.                 }
  62.             }
  63.             return *ans;
  64.         }
  65.         Matrix& operator* (double a) {
  66.             Matrix *ans = new Matrix(n, m);
  67.             for (int i = 0; i < n; i++) {
  68.                 for (int j = 0; j < m; j++) {
  69.                     ans->data[i][j] += data[i][j]*a;
  70.                 }
  71.             }
  72.             return *ans;
  73.         }
  74.         Matrix& transpose() {
  75.             Matrix *ans = new Matrix(m, n);
  76.             for (int i = 0; i < n; i++) {
  77.                 for (int j = 0; j < m; j++) {
  78.                     ans->data[j][i] = data[i][j];
  79.                 }
  80.             }
  81.             return *ans;
  82.         }
  83.         Matrix& operator| (const Matrix& matrix) {
  84.             if (matrix.n != n) {
  85.                 throw runtime_error("Error: the dimensional problem occurred");
  86.             }
  87.             Matrix *ans = new Matrix(n, matrix.m + m);
  88.             for (int i = 0; i < n; i++) {
  89.                 for (int j = 0; j < m; j++) {
  90.                     ans->data[i][j] = data[i][j];
  91.                 }
  92.                 for (int j = m; j < matrix.m + m; j++) {
  93.                     ans->data[i][j] = matrix.data[i][j-m];
  94.                 }
  95.             }
  96.             return *ans;
  97.         }
  98.         vector<double>& operator[] (int x) {
  99.             if (x >= n) {
  100.                 throw runtime_error("Error: index out of range");
  101.             }
  102.             return data[x];
  103.         }
  104.         friend istream& operator>> (istream &is, Matrix &matrix);
  105.         friend ostream& operator<< (ostream &os, Matrix &matrix);
  106.         double get(int i, int j) {
  107.             if (i >= n || j >= m) {
  108.                 throw out_of_range("Out of bound for get");
  109.             }
  110.             return data[i][j];
  111.         }
  112.         bool is_square() {
  113.             return n == m;
  114.         }
  115.         int get_n() {
  116.             return n;
  117.         }
  118.         int get_m() {
  119.             return m;
  120.         }
  121.         void normalize_diagonally() {
  122.             for (int i = 0; i < n; i++) {
  123.                 if (data[i][i] != 0) {
  124.                     for (int j = i + 1; j < m; j++) {
  125.                         data[i][j] /= data[i][i];
  126.                     }
  127.                     data[i][i] = 1;
  128.                 }
  129.             }
  130.         }
  131.         void separate() {
  132.             // for (int i = 0; i < n; i++) {
  133.             //     for (int j = 0; j < m - 1; j++) {
  134.             //         // if (abs(data[i][j]) < 0.005) {
  135.             //         //     // // cout << fixed << abs(data[i][j]) << ' ';
  136.             //         // } else
  137.             //         // // // cout << fixed << data[i][j] << ' ';
  138.             //     }
  139.             //     // cout << '\n';
  140.             // }
  141.             // for (int j = 0; j < n; j++) {
  142.             //     if (abs(data[j][m - 1]) < 0.005) {
  143.             //         // cout << fixed << abs(data[j][m - 1]) << ' ';
  144.             //     } else
  145.             //     // cout << fixed << data[j][m - 1] << ' ';
  146.             // }
  147.             // // cout << '\n';
  148.         }
  149.         void last() {
  150.             // for (int j = 0; j < n; j++) {
  151.             //     if (abs(data[j][m - 1]) < 0.005) {
  152.             //         // cout << fixed << abs(data[j][m - 1]) << ' ';
  153.             //     } else
  154.             //     // cout << fixed << data[j][m - 1] << ' ';
  155.             // }
  156.             // // cout << '\n';
  157.         }
  158.         Matrix& inverse();
  159.         operator Vector();
  160. };
  161.  
  162. class Vector : public Matrix {
  163. public:
  164.     double normalize() {
  165.         double l = 0;
  166.         for (int i = 0; i < n; i++) {
  167.             l += data[i][0] * data[i][0];
  168.         }
  169.         return l;
  170.     }
  171.     Vector(vector<double> input) : Matrix(input.size(), 1) {
  172.         for (int i = 0 ; i < input.size(); i++) {
  173.             data[i][0] = input[i];
  174.         }
  175.     }
  176.     Vector() {
  177.         n = 0;
  178.         m = 0;
  179.     }
  180.     double operator[] (int i) {
  181.         if (i >= n) {
  182.             throw runtime_error("Error: Out of range");
  183.         }
  184.         return data[i][0];
  185.     }
  186.     friend istream& operator>> (istream &is, Vector& vector);
  187. };
  188.  
  189. Matrix::operator Vector() {
  190.     if (m != 1) {
  191.         throw runtime_error("Error in cast from Matrix to Vector");
  192.     }
  193.     vector<double> a;
  194.     for (auto el : data) {
  195.         a.push_back(el[0]);
  196.     }
  197.     return *(new Vector(a));
  198. }
  199.  
  200. class LeastSquareAproxMatrix : public Matrix {
  201. public:
  202.     LeastSquareAproxMatrix(int m, vector<double> input) : Matrix(input.size(), m + 1) {
  203.         for (int i = 0; i < input.size(); i++) {
  204.             double pow = 1;
  205.             for (int j = 0; j <= m; j++) {
  206.                 data[i][j] = pow;
  207.                 pow *= input[i];
  208.             }
  209.         }
  210.     }
  211.  
  212. };
  213.  
  214. class SquareMatrix : public Matrix {
  215. public:
  216.     SquareMatrix() {
  217.         n = 0;
  218.         m = 0;
  219.     }
  220.     SquareMatrix* upper_triangular() {
  221.         SquareMatrix* ans = new SquareMatrix(n);
  222.         for (int i = 0; i < n; i++) {
  223.             for (int j = i + 1; j < n; j++) {
  224.                 (*ans)[i][j] = data[i][j];
  225.             }
  226.         }
  227.         return ans;
  228.     }
  229.     SquareMatrix* lower_triangular() {
  230.         SquareMatrix* ans = new SquareMatrix(n);
  231.         for (int i = 0; i < n; i++) {
  232.             for (int j = 0; j < i; j++) {
  233.                 (*ans)[i][j] = data[i][j];
  234.             }
  235.         }
  236.         return ans;
  237.     }
  238.     SquareMatrix& build_apha() {
  239.         SquareMatrix* ans = new SquareMatrix(n);
  240.         for (int i = 0; i < n; i++) {
  241.             for (int j = 0; j < m; j++) {
  242.                 if (i == j) {
  243.                     (*ans)[i][j] = 0;
  244.                 } else {
  245.                     (*ans)[i][j] = -data[i][j]/data[i][i];
  246.                 }
  247.             }
  248.         }
  249.         return *ans;
  250.     }
  251.     Vector& build_beta(Vector& b) {
  252.         vector<double> a;
  253.         for (int i = 0; i < n; i++) {
  254.             a.push_back(b[i] / data[i][i]);
  255.         }
  256.         Vector* ans = new Vector(a);
  257.         return *ans;
  258.     }
  259.     bool is_diagonal_dominant() {
  260.         for (int i = 0; i < n; i++) {
  261.             double sum = 0;
  262.             for (int j = 0; j < m; j++) {
  263.                 if (i != j) {
  264.                     sum += data[i][j];
  265.                 }
  266.             }
  267.             if (data[i][i] < sum) {
  268.                 return 0;
  269.             }
  270.         }
  271.         return 1;
  272.     }
  273.     SquareMatrix(int x) : Matrix(x, x) {};
  274.     friend istream& operator>> (istream &is, SquareMatrix &matrix);
  275.     void modify_for_l() {
  276.         for (int i = 0; i < n; i++) {
  277.             for (int j = 0; j <= i; j++) {
  278.                 data[i][j] = 1;
  279.             }
  280.         }
  281.     }
  282. };
  283.  
  284. class IdentityMatrix : public SquareMatrix {
  285. public:
  286.     IdentityMatrix() {
  287.         n = 0;
  288.         m = 0;
  289.     }
  290.     IdentityMatrix(int x) : SquareMatrix(x) {
  291.         for (int i = 0; i < x; i++) {
  292.             data[i][i] = 1;
  293.         }
  294.     }
  295. };
  296. class EliminationMatrix : public IdentityMatrix {
  297. public:
  298.     EliminationMatrix(Matrix& matrix, int j, int i) : IdentityMatrix(matrix.get_n()) {
  299.         data[j - 1][i - 1] = -matrix.get(j - 1, i - 1) / matrix.get(i - 1, i - 1);
  300.     }
  301.     void negative() {
  302.         for(int i = 0; i < n; i++) {
  303.             for (int j = 0; j < m; j++) {
  304.                 data[i][j] = -data[i][j];
  305.             }
  306.         }
  307.     }
  308. };
  309.  
  310. class PermutationMatrix : public IdentityMatrix {
  311. public:
  312.     PermutationMatrix(int x, int i, int j) : IdentityMatrix(x) {
  313.         data[j - 1][j - 1] = 0;
  314.         data[i - 1][i - 1] = 0;
  315.         data[i - 1][j - 1] = 1;
  316.         data[j - 1][i - 1] = 1;
  317.     }
  318. };
  319.  
  320. istream& operator>> (istream &is, SquareMatrix& matrix) {
  321.     is >> matrix.n;
  322.     matrix.m = matrix.n;
  323.     matrix.data.assign(matrix.n, vector<double>(matrix.m));
  324.     for (int i = 0; i < matrix.n; i++) {
  325.         for (int j = 0; j < matrix.m; j++) {
  326.             is >> matrix.data[i][j];
  327.         }
  328.     }
  329.     return is;
  330. }
  331. istream& operator>> (istream &is, Vector& matrix) {
  332.     is >> matrix.n;
  333.     matrix.m = 1;
  334.     matrix.data.assign(matrix.n, vector<double>(matrix.m));
  335.     for (int i = 0; i < matrix.n; i++) {
  336.         for (int j = 0; j < matrix.m; j++) {
  337.             is >> matrix.data[i][j];
  338.         }
  339.     }
  340.     return is;
  341. }
  342. istream& operator>> (istream &is, Matrix& matrix) {
  343.     is >> matrix.n >> matrix.m;
  344.     matrix.data.assign(matrix.n, vector<double>(matrix.m));
  345.     for (int i = 0; i < matrix.n; i++) {
  346.         for (int j = 0; j < matrix.m; j++) {
  347.             is >> matrix.data[i][j];
  348.         }
  349.     }
  350.     return is;
  351. }
  352. ostream& operator<< (ostream &os, Matrix& matrix) {
  353.     for (int i = 0; i < matrix.n; i++) {
  354.         for (int j = 0; j < matrix.m; j++) {
  355.             if (abs(matrix.data[i][j]) < 0.005) {
  356.                 os << fixed << abs(matrix.data[i][j]) << ' ';
  357.             } else
  358.             os << fixed << matrix.data[i][j] << ' ';
  359.         }
  360.         os << '\n';
  361.     }
  362.     return os;
  363. }
  364.  
  365. bool gaus_elimination(Matrix& a, int &step) {
  366.     bool sign = 1;
  367.     for (int i = 0; i < a.get_n(); i++) {
  368.         double m = 0;
  369.         int m_i = 0;
  370.         for (int j = i; j < a.get_n(); j++) {
  371.             if (abs(a.get(j, i)) > abs(m))  {
  372.                 m = a.get(j, i);
  373.                 m_i = j;
  374.             }
  375.         }
  376.         if (m_i != i) {
  377.             // // cout << fixed  << "step #" << step << ": permutation" << '\n';
  378.             PermutationMatrix p(a.get_n(), m_i + 1, i + 1);
  379.             a = p * a;
  380.             // a.separate();
  381.             sign = !sign;
  382.             step++;
  383.         }
  384.         for (int j = i + 1; j < a.get_n(); j++) {
  385.             if (a.get(j, i) != 0) {
  386.                 // // cout << fixed  << "step #" << step << ": elimination" << '\n';
  387.                 step++;
  388.                 EliminationMatrix e(a, j + 1, i + 1);
  389.                 a = e * a;
  390.                 // a.separate();
  391.             }  
  392.         }
  393.     }
  394.     return sign;
  395. }
  396.  
  397. void gaus_back_elimination(Matrix& a, int &step) {
  398.     for (int i = a.get_n() - 1; i >= 0; i--) {
  399.         if (a.get(i, i) != 0) {
  400.             for (int j = i - 1 ; j >= 0; j--) {
  401.                 if (a.get(j, i) != 0) {
  402.                     // // cout << fixed  << "step #" << step << ": elimination" << '\n';
  403.                     step++;
  404.                     EliminationMatrix e(a, j + 1, i + 1);
  405.                     a = e * a;
  406.                     a.separate();
  407.                 }
  408.             }
  409.         }
  410.     }
  411. }
  412.  
  413. Matrix& Matrix::inverse() {
  414.     Matrix* ans = new Matrix(n, m);
  415.     IdentityMatrix I(n);
  416.     Matrix aug = *this | I;
  417.     int step = 1;
  418.     gaus_elimination(aug, step);
  419.     gaus_back_elimination(aug, step);
  420.     aug.normalize_diagonally();
  421.     for (int j = 0; j < aug.get_n(); j++) {
  422.         for (int i = aug.get_n(); i < aug.get_m(); i++) {
  423.             (*ans)[j][i - aug.get_n()] = aug.get(j, i);
  424.         }
  425.     }
  426.  
  427.     return *ans;
  428. }
  429.  
  430. class ColumnVector : public Matrix {
  431. public:
  432.     ColumnVector() {
  433.         n = 0;
  434.         m = 0;
  435.     }
  436.     ColumnVector(int x) {
  437.         n = x;
  438.         m = 1;
  439.         data.resize(n, vector<double>(1));
  440.     }
  441.     friend istream& operator>> (istream &is, ColumnVector &matrix);
  442.     friend ostream& operator<< (ostream &os, ColumnVector &matrix);
  443. };
  444.  
  445. istream& operator>> (istream &is, ColumnVector& matrix) {
  446.     is >> matrix.n;
  447.     matrix.m = 1;
  448.     matrix.data.resize(matrix.n,vector<double>(1));
  449.     for (int i = 0; i < matrix.n; i++) {
  450.         is >> matrix.data[i][0];
  451.     }
  452.     return is;
  453. }
  454. ostream& operator<< (ostream &os, ColumnVector& matrix) {
  455.     for (int i = 0; i < matrix.n; i++) {
  456.         os << matrix.data[i][0] << ' ';
  457.     }
  458.     os << '\n';
  459.     return os;
  460. }
  461.  
  462. bool is_singular( Matrix a) {
  463.     for (int i = 0; i < a.get_n(); i++) {
  464.         double m = 0;
  465.         int m_i = 0;
  466.         for (int j = i; j < a.get_n(); j++) {
  467.             if (abs(a.get(j, i)) > abs(m))  {
  468.                 m = a.get(j, i);
  469.                 m_i = j;
  470.             }
  471.         }
  472.         if (m_i != i) {
  473.             PermutationMatrix p(a.get_n(), m_i + 1, i + 1);
  474.             a = p * a;
  475.         }
  476.         for (int j = i + 1; j < a.get_n(); j++) {
  477.             if (a.get(j, i) != 0) {
  478.                 EliminationMatrix e(a, j + 1, i + 1);
  479.                 a = e * a;
  480.             }  
  481.         }
  482.     }
  483.     for (int i = 0; i < a.get_n(); i++) {
  484.         if (a.get(i, i) == 0) {
  485.             return 0;
  486.         }
  487.     }
  488.     return 1;
  489. }
  490.  
  491. double calc_epsilon(Vector& a, Vector& b) {
  492.     if (a.get_n() != b.get_n()) {
  493.         throw runtime_error("cannot calculate epsilon");
  494.     }
  495.     double ans = 0;
  496.     for (int i = 0; i < a.get_n(); i++) {
  497.         ans += (a[i] - b[i]) * (a[i] - b[i]);
  498.     }
  499.     return sqrt(ans);
  500. }
  501.  
  502. void jacobi_method(SquareMatrix& alpha, Vector& beta, Vector& x, double epsilon, int step) {
  503.     cout << "x(" << step << "):\n";
  504.     Vector new_x = alpha*x + beta;
  505.     cout << new_x;
  506.     cout << "e:\n" << calc_epsilon(new_x, x) << '\n';
  507.     if (calc_epsilon(new_x, x) < epsilon) {
  508.         x = new_x;
  509.         return;
  510.     }
  511.     jacobi_method(alpha, beta, new_x, epsilon, step + 1);
  512.     x = new_x;
  513. }
  514.  
  515. double min(vector<double>& a) {
  516.     double m = 1000000000;
  517.     for (double el : a) {
  518.         m = min(el, m);
  519.     }
  520.     return m;
  521. }
  522. double max(vector<double>& a) {
  523.     double m = -100000000;
  524.     for (double el : a) {
  525.         m = max(el, m);
  526.     }
  527.     return m;
  528. }
  529.  
  530. void seidel_method(Matrix& B, Matrix& C, Vector& beta, Vector& x, double epsilon, int step) {
  531.     cout << "x(" << step << "):\n";
  532.     IdentityMatrix I(B.get_n());
  533.     Vector new_x = (I - B).inverse() * C * x + (I - B).inverse() * beta;
  534.     cout << new_x;
  535.     cout << "e:\n" << calc_epsilon(new_x, x) << '\n';
  536.     if (calc_epsilon(new_x, x) < epsilon) {
  537.         x = new_x;
  538.         return;
  539.     }
  540.     seidel_method(B, C, beta, new_x, epsilon, step + 1);
  541.     x = new_x;
  542. }
  543. FILE* pipe;
  544. int main() {
  545.     pipe = popen("/usr/bin/gnuplot -persist", "w");
  546.     cout.precision(2);
  547.     double v, k;
  548.     cin >> v;
  549.     cin >> k;
  550.     double a1, b1, a2, b2;
  551.     cin >> a1 >> b1 >> a2 >> b2;
  552.     double T;
  553.     cin >> T;
  554.     double n;
  555.     cin >> n;
  556.    
  557.     SquareMatrix A(2);
  558.     A[0][1] = -a2*b1/b2;
  559.     A[1][0] = a1*b2/b1;
  560.     vector<double> t;
  561.     vector<double> ks;
  562.     vector<double> vs;
  563.     Vector X_0({v, k});
  564.     // real solution
  565.     // for (double i = 0; i <= T; i += T / n) {
  566.     //     t.push_back(i);
  567.     //     complex<double> l1 = sqrt(abs(A[0][1] * A[1][0]))*1i;
  568.     //     complex<double> l2 = -sqrt(abs(A[0][1] * A[1][0]))*1i;
  569.     //     complex<double> y = (pow(M_E, l2*i) - pow(M_E, l1*i)) / (l2 - l1);
  570.     //     complex<double> x = pow(M_E, l1*i)-l1*y;
  571.     //     IdentityMatrix I(2);
  572.     //     Matrix V = I*x.real() + A*y.real();
  573.     //     Matrix Ans = V * X_0;
  574.     //     vs.push_back(Ans[0][0]);
  575.     //     ks.push_back(Ans[1][0]);
  576.     // }
  577.    
  578.     // just formula
  579.    
  580.     for (double i = 0; i <= T; i += T / n) {
  581.         vs.push_back(v*cos(sqrt(a1 * a2) * i) - k*(sqrt(a2)*b1/b2*sqrt(a1))*sin(sqrt(a1*a2)*i));
  582.         ks.push_back(v*(sqrt(a1)*b2/b1*sqrt(a2))*sin(sqrt(a1*a2)*i)+k*cos(sqrt(a1*a2)*i));
  583.         t.push_back(i);
  584.     }
  585.     cout << "t:\n";
  586.     for (double el : t) {
  587.         cout << fixed << el << ' ';
  588.     }
  589.     fputs(("plot ["+ to_string(0 - 1) + ":" + to_string(int (T + 1)) + "]["+ to_string(int(min(min(ks), min(vs)) - 1)) + ":" + to_string(int(max(max(ks), max(vs))+ 1)) + "] '-' using 1:2 title 'v(t)' with lines , '-' title 'k(t)' with lines").c_str(), pipe);
  590.     fprintf(pipe, "\n");
  591.     cout << "\nv:\n";
  592.     for (int i = 0; i < vs.size(); i++) {
  593.         cout << fixed << vs[i] << ' ';
  594.         fprintf(pipe, "%f\t%f\n", t[i], vs[i]);
  595.     }
  596.     fprintf(pipe, "%s\n", "e");
  597.     cout << "\nk:\n";
  598.     for (int i = 0; i < vs.size(); i++) {
  599.         cout << fixed << ks[i] << ' ';
  600.         fprintf(pipe, "%f\t%f\n", t[i], ks[i]);
  601.     }
  602.     fprintf(pipe, "%s\n", "e");
  603.     fputs(("plot ["+ to_string(int(min(ks) - 1)) + ":" + to_string(int(max(ks) + 1)) + "]["+ to_string(int(min(vs) - 1)) + ":" + to_string(int(max(vs) + 1)) + "] '-' using 1:2 title 'v(k)' with lines").c_str(), pipe);
  604.     fprintf(pipe, "\n");
  605.     for (int i = 0; i < vs.size(); i++) {
  606.         cout << fixed << ks[i] << ' ';
  607.         fprintf(pipe, "%f\t%f\n", ks[i], vs[i]);
  608.     }
  609.     fprintf(pipe, "%s\n", "e");
  610.     pclose(pipe);
  611. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!