SHARE
TWEET

Untitled

a guest Oct 21st, 2019 84 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <cmath>
  3. #include <complex>
  4. #include <vector>
  5. #include <fstream>
  6. using namespace std;
  7. ofstream out(“output.txt”);
  8.  
  9. complex<double> a0(0, 1), a1(6, 0), a2(3, 0), a3(20, 0), a4(10, 0); //коэффициенты исходного полинома
  10. vector<complex<double>> P = { a4, a3, a2, a1, a0 };
  11. complex<double> li(0, 1); //комплексная единица
  12.  
  13. complex<double> Horner(complex<double> z, vector<complex<double>> poly){ //функция, раскладывающая многочлен по схеме Горнера
  14.     complex<double> answer(0, 0);
  15.     for (int i = 0; i < poly.size(); i++) { answer = answer * z + poly[i]; }
  16.     return answer;
  17. }
  18.  
  19. double func(complex<double> x, vector<complex<double>> poly){ // вычисляет значение функции
  20.     complex<double> h = Horner(x, poly);
  21.     return abs(h)*abs(h);
  22. }
  23.  
  24. complex<double> CoordDescent(complex<double> z, vector<complex<double>> poly, double eps){ // координатный спуск
  25.     double a, b = 0.5; // шаг по вещественной и мнимой части
  26.     double current; //служебная переменная для вычисления функции
  27.     double F = func(z, poly);
  28.     while (F > eps){
  29.         if (F > (current = func(z + b, poly))) { z += b;  F = current; }
  30.         else if (F > (current = func(z - b, poly))) { z -= b; F = current; }
  31.         else { b /= 2; }
  32.         if (F <= eps) break;
  33.         if (F > (current = func(z + li * a, poly))) { z += li * a; F = current; }
  34.         else if (F > (current = func(z - li * a, poly))) { z -= li * a; F = current; }
  35.         else { a /= 2; }
  36.     }
  37.     return z;
  38. }
  39.  
  40. complex<double> deriv(complex<double> z, vector<complex<double>> poly){
  41.     complex<double> answer(0, 0);
  42.     double n = poly.size() - 1;
  43.     for (int i = 0; i < n; i++) { answer = answer * z + poly[i] * (n - i); }
  44.     return answer;
  45. }
  46.  
  47. double fi;
  48.  
  49. complex<double> grad(complex<double> z, vector<complex<double>> poly) { //вычисление градиента
  50.     complex<double> h = Horner(z, poly);
  51.     complex<double> dh = deriv(z, poly);
  52.     complex<double> gr = 2.0*(real(dh*conj(h)) - li * imag(dh*conj(h)));
  53.     fi = real(abs(h)*abs(h));
  54.     return gr;
  55. }
  56.  
  57. complex<double> CrushingStep(complex<double> z, vector<complex<double>> poly, double eps){ //градиентный метод с дроблением шага
  58.     double a, delta = 0.5;
  59.     double current;
  60.     complex<double> zk;
  61.     complex<double> zi = z;
  62.     complex<double> gr = grad(z, poly);
  63.     while (fi > eps){
  64.         while (true){
  65.             zk = zi - a * gr;
  66.             if (((current = func(zk, poly)) - fi) <= -a * delta * (real(gr)*real(gr) + imag(gr)*imag(gr))){
  67.                 zi = zk; fi = current; break;
  68.             }
  69.             else {
  70.                 a /= 2;
  71.             }
  72.         }
  73.         gr = grad(zi, poly);
  74.         a *= 2;
  75.     }
  76.     return zi;
  77. }
  78.  
  79. complex<double> ConstantStep(complex<double> z, vector<complex<double>> poly, double eps){ //градиентный метод с постоянным шагом
  80.     complex<double> zi = z;
  81.     complex<double> gr = grad(z, poly);
  82.     double a = 0.00005;
  83.     while (fi > eps){
  84.         zi -= a * gr;
  85.         gr = grad(zi, poly);
  86.     }
  87.     return zi;
  88. }
  89.  
  90. complex<double> GivenStep(complex<double> z, vector<complex<double>> poly, double eps){ //градиентный метод с заданным шагом
  91.     complex<double> zi = z;
  92.     complex<double> gr = grad(zi, poly);
  93.     double k = 1001;
  94.     while (fi > eps){
  95.         k++;
  96.         zi -= 1 / k * gr;
  97.         gr = grad(zi, poly);
  98.     }
  99.     return zi;
  100. }
  101.  
  102. vector<complex<double>> polynom = P;
  103.  
  104. void division(complex<double> z){ //деление полинома на корень
  105.     vector<complex<double>> temp = polynom;
  106.     polynom.clear();
  107.     polynom.push_back(temp[0]);
  108.     for (int i = 0; i < temp.size() - 2; i++) { polynom.push_back(z*polynom[i] + temp[i + 1]); }
  109. }
  110.  
  111. void SearchRoot(complex<double> x0, vector<complex<double>> poly, double eps, int method){// поиск корней
  112.     complex<double> root, n_root;
  113.     double E = 2 * eps;
  114.     auto alch = GivenStep;
  115.     polynom = P;
  116.  
  117.     switch (method){
  118.         case 1: alch = CoordDescent;  out << “Координатный спуск : “ << endl; break;
  119.         case 2: alch = CrushingStep; out << “Дробление шага : “ << endk; break;
  120.         case 3: alch = ConstantStep; out << “Постоянный шаг : “ << endl; break;
  121.         case 4: alch = GivenStep; out << “Заданный шаг : “ << endl; break;
  122.         default: break;
  123.     }
  124.  
  125.     fi = 1;
  126.     root = alch(x0, poly, eps);
  127.     division(root);
  128.  
  129.     for (int i = 1; i < poly.size() - 2; i++){
  130.         n_root = alch(x0, polynom, E);
  131.         root = alch(n_root, poly, E);
  132.         division(root);
  133.         E *= 2;
  134.     }
  135.     n_root = -polynom[1] / polynom[0];
  136.     root = alch(n_root, P, eps);
  137. }
  138.  
  139. int main(){
  140.     complex<double> point(0, 0); //точка, из которой мы вызываем методы
  141.     SearchRoot(point, P, 1e-6, 1);
  142.     SearchRoot(point, P, 1e-6, 2);
  143.     SearchRoot(point, P, 1e-6, 3);
  144.     SearchRoot(point, P, 1e-6, 4);
  145.     out.close();
  146.     return 0;
  147. }
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