SHARE
TWEET

dfp

a guest Nov 9th, 2019 80 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /*  Function Rosenbrock
  2. F(x,y)=(1-x)^2+100(y-x^2)^2
  3. partial derivative x : f = -2+2*x-400*y*x+400*x^3
  4. partial derivative y : f = 200*(y-x^2)
  5.  
  6. */
  7.  
  8. #include <conio.h>
  9. #include <math.h>
  10. #include <iostream>
  11.  
  12.  
  13.  
  14. using namespace std;
  15.  
  16.  
  17.  
  18. double ** zeroing(double **A) {
  19.     for (int i = 0; i < 2; i++) {
  20.         for (int j = 0; j < 2; j++)
  21.             A[i][j] = 0;
  22.     }
  23.     return A;
  24. }
  25.  
  26. double** product(double **A, double **U) {
  27.  
  28.     double **c = new double*[2];
  29.     for (int i = 0; i < 2; i++)
  30.         c[i] = new double[2];
  31.  
  32.     for (int i = 0; i < 2; i++) {
  33.         for (int j = 0; j < 2; j++) {
  34.             c[i][j] = 0;
  35.             for (int t = 0; t < 2; t++)
  36.                 c[i][j] += A[i][t] * U[t][j];
  37.         }
  38.     }
  39.     return c;
  40.  
  41. }
  42.  
  43. double* product_vector(double **A, double *b) {
  44.  
  45.     double *a = new double[2];
  46.  
  47.     for (int i = 0; i < 2; i++) {
  48.         a[i] = 0;
  49.         for (int j = 0; j < 2; j++)
  50.             a[i] += A[i][j] * b[j];
  51.     }
  52.     return a;
  53. }
  54.  
  55. double  deviation(double **A) {
  56.     double sum = 0, div = 0;
  57.  
  58.     for (int i = 0; i < 2; ++i) {
  59.         for (int j = 0; j < 2; ++j)
  60.             sum += A[i][j] * A[i][j];
  61.         div += sum;
  62.         sum = 0;
  63.     }
  64.  
  65.     return sqrt(div);
  66. }
  67.  
  68. int main() {
  69.  
  70.     double k=0, k_max, eps = 0.00001, t, number1, number2;
  71.     cout << "enter max approximation: ";
  72.     cin >> k_max;
  73.     cout << "enter point of min 0<t<1: ";
  74.     cin >> t;
  75.  
  76.     double **A = new double*[2]; //identity matrix
  77.     for (int i = 0; i < 2; i++)
  78.         A[i] = new double[2];
  79.  
  80.     double **B = new double*[2];
  81.     for (int i = 0; i < 2; i++)
  82.         B[i] = new double[2];
  83.  
  84.     zeroing(B);
  85.  
  86.     double **C = new double*[2];
  87.     for (int i = 0; i < 2; i++)
  88.         C[i] = new double[2];
  89.     zeroing(C);
  90.  
  91.     double **matrix1 = new double*[2];
  92.     for (int i = 0; i < 2; i++)
  93.         matrix1[i] = new double[2];
  94.  
  95.     zeroing(matrix1);
  96.  
  97.  
  98.     double **matrix2 = new double*[2];
  99.     for (int i = 0; i < 2; i++)
  100.         matrix2[i] = new double[2];
  101.  
  102.     zeroing(matrix2);
  103.  
  104.     double *a = new double[2];
  105.     for (int i = 0; i < 2; i++)
  106.         a[i] = 0;
  107.  
  108.     double *b = new double[2];
  109.     for (int i = 0; i < 2; i++)
  110.         b[i] = 0;
  111.  
  112.     double *direction = new double[2];
  113.     for (int i = 0; i < 2; i++)
  114.         direction[i] = 0;
  115.  
  116.     for (int i = 0; i < 2; i++) { // enter identity matrix
  117.         for (int j = 0; j < 2; j++) {
  118.             if (i == j)
  119.                 A[i][j] = 1;
  120.             else
  121.                 A[i][j] = 0;
  122.         }
  123.     }
  124.     double vector[2] = { 0,0 };
  125.     cout << "enter points x and y: ";
  126.     cin >> vector[0] >> vector[1];
  127.     double gradient[2] = { 0,0 };
  128.     double alpha[2] = { 0,0 };
  129.     double beta[2] = { 0,0 };
  130.     double beta_demo[2] = { 0,0 };
  131.  
  132.  
  133.     gradient[0] = -2 + 2 * vector[0] - 400 * vector[1] * vector[0] + 400 * pow(vector[0], 3);    // function's gradient evaluation at a point X
  134.     gradient[1] = 200 * (vector[1] - pow(vector[0], 2));                                        // function's gradient evaluation at a point Y
  135.     //cout << gradient[0] << " "<<gradient[1]<<endl;
  136.    
  137.     do {
  138.  
  139.         number1 = 0, number2 = 0;
  140.  
  141.         direction = product_vector(A, gradient);
  142.  
  143.         /*for (int i = 0; i < 2; i++)
  144.             cout << direction[i]<<" ";
  145.         cout << endl;*/
  146.  
  147.         direction[0] = -direction[0];
  148.         direction[1] = -direction[1];
  149.        
  150.     /*  for (int i = 0; i < 2; i++)
  151.             cout << direction[i]<<" ";
  152.             cout << endl;*/
  153.        
  154.         for (int i = 0; i < 2; i++)                        // old points
  155.             alpha[i] = vector[i];
  156.  
  157.  
  158.         /*for (int i = 0; i < 2; i++)
  159.             cout << alpha[i]<<" ";
  160.         cout << endl;
  161. */
  162.  
  163.         for (int i = 0; i < 2; i++)                        // new points
  164.             vector[i] = vector[i] + t * direction[i];
  165.  
  166.         /*for (int i = 0; i < 2; i++)
  167.             cout << vector[i] << " ";
  168.         cout << endl;*/
  169.  
  170.         for (int i = 0; i < 2; i++)                        // new points - old points
  171.             alpha[i] = vector[i] - alpha[i];
  172.  
  173.     /*  for (int i = 0; i < 2; i++)
  174.             cout << alpha[i] << " ";
  175.         cout << endl;*/
  176.  
  177.  
  178.         for (int i = 0; i < 2; i++)                       //old values of grad
  179.             beta[i] = gradient[i]; 
  180.  
  181.         /*for (int i = 0; i < 2; i++)
  182.             cout << beta[i] << " ";
  183.         cout << endl;
  184. */
  185.  
  186.         gradient[0] = -2 + 2 * vector[0] - 400 * vector[1] * vector[0] + 400 * pow(vector[0], 3);  //new values of grad
  187.         gradient[1] = 200 * (vector[1] - pow(vector[0], 2));
  188.  
  189.         /*cout << gradient[0] << " " << gradient[1] << endl;*/
  190.  
  191.         for (int i = 0; i < 2; i++)                         //new grad - old grad
  192.             beta[i] = gradient[i] - beta[i];
  193.  
  194.         //for (int i = 0; i < 2; i++)
  195.         //  cout << beta[i] << " ";
  196.         //cout << endl;
  197.  
  198.  
  199.  
  200.         for (int i = 0; i < 2; i++) {                        //begin formula
  201.             for (int j = 0; j < 2; j++) {
  202.                 matrix1[i][j] = alpha[i] * alpha[j];
  203.             }
  204.             number1 += alpha[i] * beta[i];
  205.         }
  206.         /*cout << "vot"<<endl;
  207.         for (int i = 0; i < 2; i++) {
  208.             for (int j = 0; j < 2; j++)
  209.                 cout << matrix1[i][j] << " ";
  210.             cout << endl;
  211.         }
  212.         cout << number1;*/
  213.  
  214.         for (int i = 0; i < 2; i++) {          
  215.             for (int j = 0; j < 2; j++)
  216.                 B[i][j] = matrix1[i][j] / number1;
  217.         }
  218.  
  219.        
  220.         a = product_vector(A, beta);
  221.  
  222.         for (int i = 0; i < 2; i++) {                      
  223.             for (int j = 0; j < 2; j++) {
  224.                 matrix2[i][j] = a[i] * beta[j];
  225.             }
  226.         }
  227.  
  228.         matrix2 = product(matrix2, A);
  229.  
  230.         b = product_vector(A, beta);
  231.  
  232.         for (int i = 0; i < 2; i++) {
  233.             number2 += b[i] * beta[i];
  234.         }
  235.  
  236.         for (int i = 0; i < 2; i++) {
  237.             for (int j = 0; j < 2; j++)
  238.                 C[i][j] = matrix2[i][j] / number2;
  239.         }
  240.  
  241.         for (int i = 0; i < 2; i++) {
  242.             for (int j = 0; j < 2; j++)
  243.                 A[i][j]  = A[i][j]+B[i][j]-C[i][j];           //Ak=Ak-1 + B - C;
  244.         }
  245.         k++;
  246.         cout <<endl<<"grad"<< sqrt(pow(gradient[0], 2) + pow(gradient[1], 2))<< " " << k<<endl;
  247.     } while ((sqrt(pow(gradient[0],2)+pow(gradient[1],2)) > eps) && (k < k_max));
  248.  
  249.  
  250.  
  251.     //cout << sqrt(pow(gradient[0], 2) + pow(gradient[1], 2))<< " " << k<<endl;
  252.  
  253.  
  254.  
  255.     cout << vector[0] << " " << vector[1]<<endl;
  256.     cout << pow((1 - vector[0]), 2) + 100 * pow((vector[1] - pow(vector[0], 2)), 2);
  257.  
  258.  
  259.  
  260.     _getch();
  261.     return 0;
  262. }
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