Advertisement
Guest User

Untitled

a guest
Oct 21st, 2019
116
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 4.84 KB | None | 0 0
  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. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement