SHARE
TWEET

Untitled

a guest Feb 27th, 2020 87 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <math.h>
  3. #include <vector>
  4. #include <cmath>
  5. #include <iomanip>
  6.  
  7.  
  8. void inversion(double **A, int N)
  9. {
  10.     double temp;
  11.  
  12.     double **E = new double *[N];
  13.  
  14.     for (int i = 0; i < N; i++)
  15.         E[i] = new double [N];
  16.  
  17.     for (int i = 0; i < N; i++)
  18.         for (int j = 0; j < N; j++)
  19.         {
  20.             E[i][j] = 0.0;
  21.  
  22.             if (i == j)
  23.                 E[i][j] = 1.0;
  24.         }
  25.  
  26.     for (int k = 0; k < N; k++)
  27.     {
  28.         temp = A[k][k];
  29.  
  30.         for (int j = 0; j < N; j++)
  31.         {
  32.             A[k][j] /= temp;
  33.             E[k][j] /= temp;
  34.         }
  35.  
  36.         for (int i = k + 1; i < N; i++)
  37.         {
  38.             temp = A[i][k];
  39.  
  40.             for (int j = 0; j < N; j++)
  41.             {
  42.                 A[i][j] -= A[k][j] * temp;
  43.                 E[i][j] -= E[k][j] * temp;
  44.             }
  45.         }
  46.     }
  47.  
  48.     for (int k = N - 1; k > 0; k--)
  49.     {
  50.         for (int i = k - 1; i >= 0; i--)
  51.         {
  52.             temp = A[i][k];
  53.  
  54.             for (int j = 0; j < N; j++)
  55.             {
  56.                 A[i][j] -= A[k][j] * temp;
  57.                 E[i][j] -= E[k][j] * temp;
  58.             }
  59.         }
  60.     }
  61.  
  62.     for (int i = 0; i < N; i++)
  63.         for (int j = 0; j < N; j++)
  64.             A[i][j] = E[i][j];
  65.  
  66.     for (int i = 0; i < N; i++)
  67.         delete [] E[i];
  68.  
  69.     delete [] E;
  70. }
  71.  
  72. int main() {
  73.     int n = 0;
  74.  
  75.     std::cout << "Enter matrix dimension (n):" << std::endl;
  76.     std::cin >> n;
  77.  
  78.     float *diag_a = new float[n];
  79.     float *diag_b = new float[n];
  80.     float *diag_c = new float[n];
  81.     float *vect_d = new float[n];
  82.     float *vect_x = new float[n];
  83.  
  84.     float *vect_d1 = new float[n];
  85.     float *vect_r = new float[n];
  86.     float *vect_e = new float[n];
  87.  
  88.     float *alpha = new float[n];
  89.     float *beta = new float[n];
  90.  
  91. //    int N;
  92. //    std::cout << "Enter N: ";
  93. //    std::cin >> N;
  94.     double  **A = new double *[n];
  95.  
  96.     std::cout << "Enter main diagonal (a):" << std::endl;
  97.     for (int i = 0; i < n; i++) {
  98.         std::cin >> diag_a[i];
  99.     }
  100.  
  101.     std::cout << "Enter side diagonal 1 (b):" << std::endl;
  102.     for (int i = 0; i < n - 1; i++) {
  103.         std::cin >> diag_b[i];
  104.     }
  105.  
  106.     std::cout << "Enter side diagonal 2 (c):" << std::endl;
  107.     for (int i = 0; i < n - 1; i++) {
  108.         std::cin >> diag_c[i];
  109.     }
  110.  
  111.     std::cout << "Enter side free vector (d):" << std::endl;
  112.     for (int i = 0; i < n; i++) {
  113.         std::cin >> vect_d[i];
  114.     }
  115.  
  116.  
  117.     //Check diagonal prevalence conditions
  118.  
  119.     for (int i = 1; i < n - 1; i++) {
  120.         if (abs(diag_a[i]) < abs(diag_c[i - 1] + diag_b[i])) {
  121.             std::cout << "The first condition of diagonal dominance is not satisfied!" << std::endl;
  122.             return 0;
  123.         } else std::cout << "Ok. The first condition of diagonal dominance is satisfied!" << std::endl;
  124.  
  125.  
  126.         if ((abs(diag_b[i]) / abs(diag_a[i])) > 1) {
  127.             std::cout << "The second condition of diagonal dominance is not satisfied!" << std::endl;
  128.             return 0;
  129.         } else std::cout << "Ok. The second condition of diagonal dominance is satisfied!" << std::endl;
  130.  
  131.         if ((abs(diag_c[i - 1]) / abs(diag_b[i])) > abs(diag_c[i - 1] + diag_b[i])) {
  132.             std::cout << "The third condition of diagonal dominance is not satisfied!" << std::endl;
  133.             return 0;
  134.         } else std::cout << "Ok. The third condition of diagonal dominance is satisfied!" << std::endl;
  135.     }
  136.  
  137.  
  138.     //Find vector x
  139.     if (diag_a[0] != 0) {
  140.         alpha[0] = -(diag_b[0] / diag_a[0]);
  141.         beta[0] = vect_d[0] / diag_a[0];
  142.     } else return 0;
  143.  
  144.  
  145.     for (int i = 1; i < n; i++) {
  146.         alpha[i] = -diag_b[i] / (alpha[i - 1] * diag_c[i - 1] + diag_a[i]);
  147.         beta[i] = (vect_d[i] - diag_c[i - 1] * beta[i - 1]) / (alpha[i - 1] * diag_c[i - 1] + diag_a[i]);
  148.     }
  149.  
  150.     vect_x[n - 1] = beta[n - 1];
  151.  
  152.     for (int i = n - 2; i >= 0; i--) {
  153.         vect_x[i] = (alpha[i] * vect_x[i + 1] + beta[i]);
  154.         //vect_x[i] = (vect_d[i] - diag_c[i - 1] * beta[i - 1]) / (diag_a[i] + (diag_c[i - 1] * alpha[i - 1]));
  155.     }
  156.  
  157.     std::cout << std::endl;
  158.     std::cout << "Vector x (x*)" << std::endl;
  159.  
  160.     for (int i = 0; i < n; i++) {
  161.         std::cout << std::fixed << std::setprecision(15) << vect_x[i] << std::endl;
  162.     }
  163.  
  164. //    std::cout << std::endl;
  165. //    std::cout << "Vector d*" << std::endl;
  166. //
  167. //    for (int i = 0; i < n; i++) {
  168. //        vect_d1[i] = diag_a[i] * vect_x[i];
  169. //        std::cout << std::fixed << std::setprecision(15) << vect_d1[i] << std::endl;
  170. //
  171. //    }
  172.  
  173.  
  174.  
  175.     std::cout << "Enter matrix A:" << std::endl;
  176.  
  177.     for (int i = 0; i < n; i++) {
  178.         A[i] = new double[n];
  179.     }
  180.  
  181.     for (int i = 0; i < n; i++) {
  182.         for (int j = 0; j < n; j++) {
  183.             std::cin >> A[i][j];
  184.         }
  185.     }
  186.  
  187.     std::cout << std::endl;
  188.     for (int i = 0; i < n; i++) {
  189.         for (int j = 0; j < n; j++) {
  190.             vect_d1[i] = vect_d1[i] + (A[i][j] * vect_x[i]);
  191.         }
  192.         //std::cout << vect_d1[i] << std::endl;
  193.         std::cout << std::fixed << std::setprecision(15) << "d*" << vect_d1[i] << std::endl;
  194.  
  195.     }
  196.  
  197.     std::cout << std::endl;
  198.     std::cout << "Vector r" << std::endl;
  199.  
  200.     for (int i = 0; i < n; i++) {
  201.         vect_r[i] = vect_d[i] - vect_d1[i];
  202.         std::cout << std::fixed << std::setprecision(15) << vect_r[i] << std::endl;
  203.  
  204.     }
  205.  
  206.     std::cout << "Matrix A^-1:" << std::endl;
  207.     inversion(A, n);
  208.     for (int i = 0; i < n; i++)
  209.     {
  210.         for (int j = 0; j < n; j++)
  211.             std::cout << A[i][j] << "  ";
  212.  
  213.         std::cout << std::endl;
  214.     }
  215.  
  216.     std::cout << std::endl;
  217.     std::cout << "Vector e" << std::endl;
  218.  
  219.     std::cout << std::endl;
  220.     for (int i = 0; i < n; i++) {
  221.         for (int j = 0; j < n; j++) {
  222.             vect_e[i] = vect_e[i] + (A[i][j] * vect_r[i]);
  223.         }
  224.         //std::cout << vect_d1[i] << std::endl;
  225.         std::cout << std::fixed << std::setprecision(15) << "e" << vect_e[i] << std::endl;
  226.  
  227.     }
  228. }
RAW Paste Data
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
Top