constk

Solving non linear equations 1.1

Sep 21st, 2020 (edited)
149
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #define _CRT_SECURE_NO_WARNINGS
  2. #include <stdio.h>
  3. #include <cmath>
  4. #include <vector>
  5. #include <stdexcept>
  6.  
  7. // Для возведения в квадрат
  8. #define sqr(a) ((a)*(a))
  9.  
  10.  
  11.  
  12. // Отрезок [a;b] на котором находится один корень уравнения
  13. struct Section {
  14.     double leftBorder; // a
  15.     double rightBorder;// b
  16. };
  17.  
  18.  
  19.  
  20. // Точность эпсилон
  21. const double EPS = 0.001;
  22.  
  23.  
  24.  
  25. // Прототипы функций для уравнений
  26. double f1(double x); // [0 ; 1] root is 0.1
  27. double f2(double x); // just y = x
  28. double f3(double x); // [1.5 ; 3.14] root is 2
  29. double f4(double x); // just y = sin(x)
  30. double f5(double x); // just x*x
  31.  
  32. // Тип указателя на функции выше
  33. typedef double (*pointFunc)(double x);
  34.  
  35.  
  36.  
  37. // Прототипы функций для вычисления производных
  38. double derivative(pointFunc f, double x);
  39. double doubleDerivative(pointFunc f, double x);
  40.  
  41. // Прототип функции sign() - возвращает знак (+1/-1)
  42. int sign(double x);
  43.  
  44.  
  45.  
  46. // Прототипы функций для решения уравнения
  47. double halfDivisionMethod(pointFunc f, const Section & section);
  48. double chordMethod(pointFunc f, const Section & section);
  49. double newtonMethod(pointFunc f, const Section & section);
  50. double interationMethod(pointFunc f, const Section & section);
  51.  
  52. // Тип указателя на функции выше
  53. typedef double(*pointMethod)(pointFunc f, const Section & section);
  54.  
  55.  
  56.  
  57. int main() {
  58.     // Здесь будут хранится отрезки, содержащие один корень уравнения
  59.     std::vector<Section> sections;
  60.  
  61.     // Здесь будет храниться найденные корни уравнения
  62.     std::vector<double> roots;
  63.  
  64.  
  65.  
  66.     //############################################################################
  67.     //# Ввод исходных данных и разбиение функции на интервалы (отделение корней) #
  68.     //############################################################################
  69.  
  70.     const pointFunc f = f1; // Здесь задаём уравнение через указатель на функцию
  71.     int n = 0;
  72.     double leftBorder = 0.0;
  73.     double rightBorder = 0.0;
  74.  
  75.     // Вводим область определения, на которой хотим искать корни
  76.     // Также вводим кол-во отрезков, необходимое для поиска всех корней
  77.     initInput:
  78.     printf("Input left border: ");
  79.     scanf("%lf", &leftBorder);
  80.     printf("Input right border: ");
  81.     scanf("%lf", &rightBorder);
  82.     printf("Input amount of sections: ");
  83.     scanf("%d", &n);
  84.     if (rightBorder <= leftBorder) {
  85.         puts("Error: incorrect input of section borders! Try again!\n");
  86.         goto initInput;
  87.     }
  88.     if (n <= 0) {
  89.         puts("Error: incorrect input amount of sections! Try again!\n");
  90.         goto initInput;
  91.     }
  92.  
  93.     // Запоминаем те отрезки [a; b], где есть корень
  94.     // Таковыми будем считать все отрезки, знак функции на концах которых различается
  95.     // Т.е. f(a) * f(b) = -1
  96.     double x = leftBorder;
  97.     double h = (rightBorder - leftBorder) / double(n);
  98.     while (x < rightBorder)
  99.     {
  100.         if (f(x) * f(x + h) <= 0)
  101.             sections.push_back({ x, x + h });
  102.         x += h;
  103.     }
  104.     if (sections.empty()) { // Если таких отрезков нет
  105.         printf("Error: something went wrong trying to get sections.");
  106.         printf("Ensure equation have roots or increase N.");
  107.         getchar();
  108.         getchar();
  109.         exit(1);
  110.     }
  111.     puts("");
  112.     //sections.push_back({ -1.0, 1.0 });
  113.  
  114.  
  115.  
  116.     //############################################
  117.     //# Выбор метода решения и вычисление корней #
  118.     //############################################
  119.  
  120.     // Запрашиваем у пользователя номер метода
  121.     pointMethod method = NULL;
  122.     backToChoice:
  123.     puts("Choose method: ");
  124.     puts("Print 1 to use half division method.");
  125.     puts("Print 2 to use chord method.");
  126.     puts("Print 3 to use Newton method.");
  127.     puts("Print 4 to use iteration method.");
  128.     int choice = 0;
  129.     scanf("%d", &choice);
  130.  
  131.     // Инициализируем указатель на функцию выбранным методом
  132.     switch (choice)
  133.     {
  134.     case 1:
  135.         method = halfDivisionMethod;
  136.         break;
  137.     case 2:
  138.         method = chordMethod;
  139.         break;
  140.     case 3:
  141.         method = newtonMethod;
  142.         break;
  143.     case 4:
  144.         method = interationMethod;
  145.         break;
  146.     default:
  147.         puts("There is no such method! Choose again!\n");
  148.         goto backToChoice;
  149.     }
  150.  
  151.     // Решаем уравнение выбранном методом на заданных отрезках
  152.     while (!sections.empty()) {
  153.         try {
  154.             // Запоминаем вычисленные корни
  155.             roots.push_back(method(f, sections.back()));
  156.         }
  157.         catch (std::runtime_error& e) { // Если что-то пошло не так
  158.             printf("Error: %s\n\n", e.what());
  159.         }
  160.         sections.pop_back();
  161.     }
  162.     puts("");
  163.  
  164.  
  165.  
  166.     //####################################
  167.     //# Вывод результатов и тестирование #
  168.     //####################################
  169.  
  170.     // Выводим все корни на экран
  171.     printf("Roots: ");
  172.     for (int i = 0; i != roots.size(); ++i)
  173.     printf("%f   ", roots.at(i));
  174.     puts("\n");
  175.  
  176.     // Проверяем корни на правильность
  177.     double currentRoot = 0.0;
  178.     double currentAnswer = 0.0;
  179.     while (!roots.empty()) {
  180.         currentRoot = roots.back();
  181.         currentAnswer = f(currentRoot);
  182.         if (sqr(currentAnswer) > EPS*2)
  183.             printf("Test failed: %f is not root of the equation! f(%f) = %f\n", currentRoot, currentRoot, f(currentRoot));
  184.         else
  185.             printf("Test passed.\n");
  186.         roots.pop_back();
  187.     }
  188.     puts("");
  189.  
  190.  
  191.  
  192.     //###################
  193.     //# Конец программы #
  194.     //###################
  195.  
  196.     puts("End of program. Press enter to exit.");
  197.     getchar();
  198.     getchar();
  199.     return 0;
  200. }
  201.  
  202.  
  203.  
  204. //############################################
  205. //# Функции f(x) для уравнений вида f(x) = 0 #
  206. //############################################
  207.  
  208. double f1(double x) {
  209.     return exp(x) - 10 * x;
  210. }
  211.  
  212. double f2(double x) {
  213.     return 2 * x;
  214. }
  215.  
  216. double f3(double x) {
  217.     return x * x - 5 * sin(x);
  218. }
  219.  
  220. double f4(double x)
  221. {
  222.     return sin(x);
  223. }
  224.  
  225. double f5(double x)
  226. {
  227.     return x*x;
  228. }
  229.  
  230.  
  231. //######################################
  232. //# Функции для нахождения производных #
  233. //######################################
  234.  
  235. // Первая производная функции в точке
  236. double derivative(pointFunc f, double x)
  237. {
  238.     //const double h = 0.0005;
  239.     //return (f(x) - f(x + h)) / h;
  240.     return exp(x) - 10; // For f1
  241.     //return 2 * x - 5 * cos(x); // For f3
  242.     //return 2 * x; // For f5
  243. }
  244.  
  245. // Вторая производная функции в точке
  246. double doubleDerivative(pointFunc f, double x)
  247. {
  248.     //const double h = 0.0005;
  249.     //return (f(x + h) - 2 * f(x) + f(x - h)) / (h * h);
  250.     return exp(x); // For f1
  251.     //return 2 + 5 * sin(x); // For f3
  252.     //return 2.0; // For f5
  253. }
  254.  
  255. // Знак (проверяет знак числа, возвращает +1 или -1)
  256. int sign(double x)
  257. {
  258.     return (x > 0) ? (1) : (-1);
  259. }
  260.  
  261. //######################################################################
  262. //# Функции для нахождения одного корня уравнения на указанном отрезке #
  263. //######################################################################
  264.  
  265. // Метод половинного деления
  266. double halfDivisionMethod(pointFunc f, const Section & section) {
  267.     int iterationsCounter = 0;
  268.     double a = section.leftBorder;
  269.     double b = section.rightBorder;
  270.     double x = (a + b) / 2.0;
  271.     double yLeft = f(a), yRight = f(b), yCenter = f(x);
  272.  
  273.     while (fabs(yRight - yLeft) > EPS  && iterationsCounter != 10000) {
  274.         if (yCenter * yRight < 0) {
  275.             a = x;
  276.             yLeft = f(a);
  277.         }
  278.         else {
  279.             b = x;
  280.             yRight = f(b);
  281.         }
  282.  
  283.         x = (a + b) / 2.0;
  284.         yCenter = f(x);
  285.         ++iterationsCounter;
  286.     }
  287.  
  288.     if (iterationsCounter == 10000)
  289.         throw std::runtime_error("Half division method does not respond.");
  290.  
  291.     return x;
  292. }
  293.  
  294. // Метод хорд
  295. double chordMethod(pointFunc f, const Section & section) {
  296.     int iterationsCounter = 0;
  297.     double a = section.leftBorder;
  298.     double b = section.rightBorder;
  299.     double result = 0.0;
  300.     double yLeft = f(a), yRight = f(b);
  301.  
  302.  
  303.     if ((yLeft*doubleDerivative(f ,a)) > 0)
  304.     {
  305.         while (fabs(b - ((b - a) / (yRight - yLeft))*yRight - b) > EPS && iterationsCounter != 100000)
  306.         {
  307.             if (sqr(yRight - yLeft) <= EPS)
  308.                 throw std::runtime_error("Chord method does not converge.");
  309.  
  310.             b = b - ((b - a) / (yRight - yLeft))*yRight;
  311.             ++iterationsCounter;
  312.             yRight = f(b);
  313.         }
  314.         result = b;
  315.     }
  316.     else
  317.     {
  318.         if (sqr(yRight - yLeft) <= EPS)
  319.             throw std::runtime_error("Chord method does not converge.");
  320.  
  321.         while (fabs(a - ((b - a) / (yRight - yLeft))*yLeft - a) > EPS && iterationsCounter != 100000)
  322.         {
  323.             a = a - ((b - a) / (yRight - yLeft))*yLeft;
  324.             ++iterationsCounter;
  325.             yLeft = f(a);
  326.         }
  327.         result = a;
  328.     }
  329.  
  330.  
  331.     if (iterationsCounter == 100000)
  332.         throw std::runtime_error("Chord method does not respond.");
  333.  
  334.     return result;
  335. }
  336.  
  337. // Метод Ньютона
  338. double newtonMethod(pointFunc f, const Section & section)
  339. {
  340.     int iterationsCounter = 0;
  341.     double a = section.leftBorder;
  342.     double b = section.rightBorder;
  343.     double yLeft = f(a), yRight = f(b);
  344.  
  345.  
  346.     //if(sign(derivative(f, a)) != sign(derivative(f, b)) || sign(doubleDerivative(f, a)) != sign(doubleDerivative(f, b)))
  347.     //  throw std::runtime_error("Newton method does not converge.");
  348.  
  349.  
  350.     double x = a;
  351.     if (!(yLeft * doubleDerivative(f, x) > 0))
  352.     {
  353.         x = b;
  354.         if (!(yRight * doubleDerivative(f, x) > 0))
  355.         {
  356.             x = (a + b) / 2;
  357.             if (!(f(x) * doubleDerivative(f, x) > 0))
  358.                 puts("Warning: lost of accuracy is possible using Newton Method!\n");
  359.         }
  360.     }
  361.  
  362.  
  363.     double next = x;
  364.     do
  365.     {
  366.         x = next;
  367.         if(sqr(derivative(f, x)) <= EPS)
  368.             throw std::runtime_error("Newton method does not converge.");
  369.  
  370.         next = x - f(x) / derivative(f, x);
  371.     } while (fabs(next - x) > EPS && iterationsCounter != 100000);
  372.  
  373.  
  374.     if (iterationsCounter == 100000)
  375.         throw std::runtime_error("Newton method does not respond.");
  376.  
  377.  
  378.     return next;
  379. }
  380.  
  381. // Метод итераций
  382. double interationMethod(pointFunc f, const Section & section)
  383. {
  384.     int iterationsCounter = 0;
  385.     double x = section.leftBorder;
  386.     double k = 0.1;
  387.  
  388.  
  389.     double next = x;
  390.     do
  391.     {
  392.         k = 0.1;
  393.         while (fabs(1 + k * derivative(f, x)) > 0.5)
  394.             k /= 10;
  395.  
  396.         x = next;
  397.         next = x + k * f(x);
  398.     } while (fabs(next - x) > EPS && iterationsCounter != 100000);
  399.  
  400.  
  401.     if (iterationsCounter == 100000)
  402.         throw std::runtime_error("Iteration method does not respond.");
  403.  
  404.  
  405.     return next;
  406. }
  407.  
RAW Paste Data