SHARE
TWEET

least-square method

frustration Oct 21st, 2019 (edited) 81 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include<iostream>
  2. #include<math.h>
  3. #include <conio.h>
  4.  
  5.  
  6. using namespace std;
  7.  
  8.  
  9. double  deviation(int k, double **A) {
  10.     double sum = 0, div=0;
  11.  
  12.     for (int i = 0; i < k + 1; ++i) {            
  13.         for (int j = 0; j < k + 1; ++j)
  14.             sum += A[i][j] * A[i][j];
  15.         div += sum;
  16.         sum = 0;
  17.     }
  18.  
  19.     return sqrt(div);
  20. }
  21.  
  22. double ** zeroing(int k, double **A) {
  23.     for (int i = 0; i < k + 1; i++) {
  24.         for (int j = 0; j < k + 1; j++)
  25.             A[i][j] = 0;
  26.     }
  27.     return A;
  28. }
  29.  
  30. double** product(int k, double **A, double **U) {
  31.    
  32.     double **c = new double*[k + 1];
  33.     for (int i = 0; i < k + 1; i++)
  34.         c[i] = new double[k + 1];
  35.  
  36.     for (int i = 0; i < k + 1; i++) {
  37.         for (int j = 0; j < k + 1; j++) {
  38.         c[i][j] = 0;
  39.             for (int t = 0; t < k + 1; t++)
  40.             c[i][j] += A[i][t] * U[t][j];
  41.         }
  42.     }
  43.     return c;
  44.  
  45. }
  46.  
  47.  
  48. int main()
  49. {
  50.     int n, k, step;
  51.     cout << "enter steps: ";
  52.     cin >> step;
  53.     cout << "enter number of points N: ";
  54.     cin >> n;
  55.     cout << "enter degree K: ";
  56.     cin >> k;
  57.     double *x = new double[n];
  58.     double *y = new double[n];
  59.     cout << "enter X: ";
  60.     for (int i = 0; i < n; ++i)
  61.         cin >> x[i];
  62.  
  63.     cout << "enter Y: ";
  64.     for (int i = 0; i < n; ++i)
  65.         cin >> y[i];
  66.  
  67.  
  68.     double **A = new double*[k + 1]; // matrix A
  69.     for (int i = 0; i < k + 1; i++)
  70.         A[i] = new double[k + 1];
  71.  
  72.     double *b = new double[k + 1]; // vector b
  73.     for (int i = 0; i < k + 1; i++)
  74.         b[i] = 0;
  75.  
  76.     double *a = new double[k + 1]; // vector a
  77.     for (int i = 0; i < k + 1; i++)
  78.         a[i] = 0;
  79.        
  80.     double **F = new double*[k + 1]; //matrix Psi
  81.     for (int i = 0; i < k + 1; i++)
  82.         F[i] = new double[k + 1];
  83.  
  84.     double **E = new double*[k + 1]; //identity matrix
  85.     for (int i = 0; i < k + 1; i++)
  86.         E[i] = new double[k + 1];
  87.    
  88.     double **c = new double*[k + 1];
  89.     for (int i = 0; i < k + 1; i++)
  90.         c[i] = new double[k + 1];
  91.  
  92.     A = zeroing(k, A);
  93.  
  94.     for (int i = 0; i < k + 1; i++) { //enter matrix A
  95.         for (int j = 0; j < k + 1; j++)
  96.             for (int t = 0; t < n; t++)
  97.                 A[i][j] += pow(x[t], i + j);       
  98.     }
  99.  
  100.  
  101.     for (int i = 0; i < k + 1; ++i) { //enter vector b
  102.         for (int j = 0; j < n; ++j) {
  103.             b[i] += pow(x[j], i)*y[j];
  104.         }
  105.     }
  106.  
  107.     double **U = new double*[k + 1]; //approximation.
  108.     for (int i = 0; i < k + 1; i++)
  109.         U[i] = new double[k + 1];
  110.  
  111.     c = zeroing(k, c);
  112.  
  113.     c = product(k, A, A);
  114.  
  115.     for (int i = 0; i < k + 1; ++i) { // first approximation
  116.         for (int j = 0; j < k + 1; ++j)
  117.             U[i][j] = A[i][j] / deviation(k,c);
  118.     }
  119.  
  120.     for (int i = 0; i < k + 1; i++) { // enter identity matrix
  121.         for (int j = 0; j < k + 1; j++) {
  122.             if (i == j)
  123.                 E[i][j] = 1;
  124.             else
  125.                 E[i][j] = 0;
  126.         }
  127.     }
  128.  
  129.     double eps = 0.00001;
  130.     int oks = 0;
  131.  
  132.  
  133.     do {                                    //Shultz's method
  134.  
  135.         c = zeroing(k, c);
  136.  
  137.         c = product(k, A, U);
  138.        
  139.         for (int i = 0; i < k + 1; ++i) {
  140.             for (int j = 0; j < k + 1; ++j)
  141.                 F[i][j] = E[i][j] - c[i][j];             //find Psi
  142.         }
  143.        
  144.         for (int i = 0; i < k + 1; ++i) {
  145.             for (int j = 0; j < k + 1; ++j) {
  146.                 c[i][j] = 0;
  147.                 c[i][j] = E[i][j] + F[i][j];
  148.             }
  149.         }
  150.  
  151.         U = product(k, U, c);
  152.  
  153.         oks++;
  154.  
  155.     } while ((oks < step ) and (deviation(k, F) > eps));
  156.  
  157.     for (int i = 0; i < k + 1; i++) {
  158.         a[i] = 0;
  159.         for (int j = 0; j < k + 1; j++)
  160.             a[i] += U[i][j] * b[j];                        //find vector a by applying inverse matrix U and vector b
  161.     }
  162.  
  163.     cout << "vector a"<<endl;
  164.     for (int i = 0; i < k + 1; i++) {                    // output a
  165.         cout << a[i] << " ";
  166.     }
  167.  
  168.     _getch();
  169.     return 0;
  170. }
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