SeriousVenom

Lab3

Apr 22nd, 2021 (edited)
603
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <cmath>
  3.  
  4. using std::cout;
  5. using std::cin;
  6. using std::endl;
  7.  
  8. static const double h = 0.1;
  9.  
  10. double derivative(double x, double y); // Дифференциальное уравнение
  11. double correct(double x_last, double y_last, double x_next, double y_cur, double h, double eps); // Проверка предсказанного значения \\ Мод. метод Эйлера
  12. double rungeK(double h, double x, double y, int number); // Вычисление k1, k2, k3, k4 \\ Метод Рунге-Кутта
  13. double yt(double x, double y);
  14. void eulerMet(bool select); // Метод Эйлера
  15. void rungeKutt(bool select); // Метод Рунге-Кутта
  16. void milnMet(); // Метод Милна
  17.  
  18. int main()
  19. {
  20.     unsigned short int pck = -1;
  21.     bool check_task = false;
  22.     cout << "Solution of the Cauchy problem" << endl;
  23.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  24.     cout << "Choose a method:\n1.Euler's method\n2.Runge-Kutt method\n3.Miln method\nEnter:";
  25.  
  26.     while (check_task == false)
  27.     {
  28.         cin >> pck;
  29.         switch (pck)
  30.         {
  31.         case 1:
  32.         {
  33.             check_task = true;
  34.             double a, b, eps;
  35.             bool mod;
  36.             cout << "0 - Not modified\n1 - Modified\nEnter: ";
  37.             cin >> mod;
  38.             if (mod)
  39.             {
  40.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  41.                 cout << "Modified Euler's method" << endl;
  42.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  43.             }
  44.             else
  45.             {
  46.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  47.                 cout << "Not modified Euler's method" << endl;
  48.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  49.             }
  50.  
  51.             eulerMet(mod);
  52.  
  53.             break;
  54.         }
  55.         case 2:
  56.         {
  57.             check_task = true;
  58.             double a, b, eps;
  59.             bool mod;
  60.             cout << "0 - Without precision\n1 - With a precision\nEnter: ";
  61.             cin >> mod;
  62.             if (mod)
  63.             {
  64.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  65.                 cout << "Runge-Kutt method with a given precision" << endl;
  66.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  67.             }
  68.             else
  69.             {
  70.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  71.                 cout << "Runge-Kutt method without precision" << endl;
  72.                 cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  73.             }
  74.  
  75.             rungeKutt(mod);
  76.  
  77.             break;
  78.         }
  79.         case 3:
  80.         {
  81.             cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  82.             cout << "\tMiln Method" << endl;
  83.             cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  84.  
  85.             milnMet();
  86.  
  87.             break;
  88.         }
  89.         default: {cout << "#Error#\nRe-enter: "; }
  90.         }
  91.     }
  92.  
  93.     return 0;
  94. }
  95.  
  96. double derivative(double x, double y)
  97. {
  98.     return (exp(x) + exp(-x));
  99. }
  100. double yt(double x)
  101. {
  102.     return(exp(x) - exp(-x));
  103. }
  104. double correct(double x_last, double y_last, double x_next, double y_cur, double h, double eps)
  105. {
  106.     double y_correct = y_cur;
  107.     do
  108.     {
  109.         y_cur = y_correct;
  110.         y_correct = y_last + h * 0.5 * (derivative(x_last, y_last) + derivative(x_next, y_cur));
  111.     } while (fabs(y_correct - y_cur) > eps);
  112.  
  113.     return y_correct;
  114. }
  115. void eulerMet(bool select)
  116. {
  117.     double a = 0, b = 1;
  118.     double n = (b - a) / h;
  119.     double* massX = new double[n + 1];
  120.     double* massY = new double[n + 1];
  121.     double* massT = new double[n + 1];
  122.     cout << "a (left border) = " << a << "\nb (right border) = " << b << endl;
  123.     cout << "h (step) = " << h << endl;
  124.     massX[0] = a;
  125.     massY[0] = 0;
  126.     switch (select)
  127.     {
  128.     case true:
  129.     {
  130.         double eps;
  131.         double x1;
  132.         double y_correct;
  133.         cout << "Enter eps = "; cin >> eps;
  134.         for (unsigned short int i = 1; i <= n; i++)
  135.         {
  136.             massX[i] = a + i * h;
  137.             massY[i] = massY[i - 1] + h * derivative(massX[i - 1], massY[i - 1]);
  138.             x1 = massX[i] + h;
  139.             massY[i] = correct(massX[i - 1], massY[i - 1], x1, massY[i], h, eps);
  140.  
  141.         }
  142.  
  143.         break;
  144.     }
  145.     case false:
  146.     {
  147.         for (unsigned short int i = 1; i <= n; i++)
  148.         {
  149.             massX[i] = a + i * h;
  150.             massY[i] = massY[i - 1] + h * derivative(massX[i - 1], massY[i - 1]);
  151.         }
  152.  
  153.         break;
  154.     }
  155.     }
  156.     for (unsigned short int i = 0; i <= n; i++)
  157.     {
  158.         massT[i] = yt(massX[i]);
  159.     }
  160.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  161.     cout << "Equation: y` = e^x + e^(-x); y(0) = 0; yt = e^x - e^(-x)" << endl;
  162.     cout << "Results: " << endl;
  163.     for (unsigned short int i = 0; i <= n; i++)
  164.     {
  165.         cout.precision(9);
  166.         cout << "X[" << i << "] = " << massX[i] << "\tY[" << i << "] = " << massY[i] << "\terror = " << fabs(massY[i] - massT[i]) << endl;
  167.     }
  168.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  169.     delete[] massX;
  170.     delete[] massY;
  171. }
  172. void rungeKutt(bool select)
  173. {
  174.     double a = 0.0, b = 1.0;
  175.     double n = (b - a) / h;
  176.     double* massX = new double[n + 1];
  177.     double* massY = new double[n + 1];
  178.     double* massT = new double[n + 1];
  179.     double k1, k2, k3, k4;
  180.     cout << "a (left border) = " << a << "\nb (right border) = " << b << endl;
  181.     cout << "h (step) = " << h << endl;
  182.     massX[0] = a;
  183.     massY[0] = 0;
  184.     switch (select)
  185.     {
  186.     case true:
  187.     {
  188.         double eps;
  189.         double h0 = (b - a) / n;
  190.         double h1 = h0 / 2.0;
  191.         double y1, y2, yv;
  192.         double sum1 = 0, sum2 = 0, sum3 = 0;
  193.         cout << "Enter eps = "; cin >> eps;
  194.         for (unsigned short int i = 1; i <= n; i++)
  195.         {
  196.             massX[i] = a + i * h0;
  197.             for (unsigned short int j = 1; i < 5; i++)
  198.             {
  199.                 sum1 += h0 * rungeK(h0, massX[i - 1], massY[i - 1], j);
  200.             }
  201.             y1 = massY[i - 1] + sum1;
  202.             for (unsigned short int j = 1; i < 5; i++)
  203.             {
  204.                 sum2 += h1 * rungeK(h1, massX[i - 1], massY[i - 1], j);
  205.             }
  206.             yv = massY[i - 1] + sum2;
  207.             for (unsigned short int j = 1; i < 5; i++)
  208.             {
  209.                 sum3 += h1 * rungeK(h1, massX[i - 1], yv, j);
  210.             }
  211.             y2 = yv + sum3;
  212.  
  213.             if (fabs(y1 - y2) < eps)
  214.             {
  215.                 for (unsigned short int i = 1; i <= n; i++)
  216.                 {
  217.                     massX[i] = a + i * h0;
  218.                     massY[i] = massY[i - 1] + (((rungeK(h0, massX[i - 1], massY[i - 1], 1)) + 2 * (rungeK(h0, massX[i - 1],
  219.                         massY[i - 1], 2) + rungeK(h0, massX[i - 1], massY[i - 1], 3)) + rungeK(h0, massX[i - 1], massY[i - 1], 4)) / 6.0);
  220.                 }
  221.             }
  222.             else
  223.             {
  224.                 for (unsigned short int i = 1; i <= n; i++)
  225.                 {
  226.                     massX[i] = a + i * h0;
  227.                     massY[i] = massY[i - 1] + (((rungeK(h1, massX[i - 1], massY[i - 1], 1)) + 2 * (rungeK(h1, massX[i - 1],
  228.                         massY[i - 1], 2) + rungeK(h1, massX[i - 1], massY[i - 1], 3)) + rungeK(h1, massX[i - 1], massY[i - 1], 4)) / 6.0);
  229.                 }
  230.             }
  231.         }
  232.         break;
  233.     }
  234.     case false:
  235.     {
  236.         for (unsigned short int i = 1; i <= n; i++)
  237.         {
  238.             massX[i] = a + i * h;
  239.             massY[i] = massY[i - 1] + (((rungeK(h, massX[i - 1], massY[i - 1], 1)) + 2 * (rungeK(h, massX[i - 1],
  240.                 massY[i - 1], 2) + rungeK(h, massX[i - 1], massY[i - 1], 3)) + rungeK(h, massX[i - 1], massY[i - 1], 4)) / 6.0);
  241.         }
  242.         break;
  243.     }
  244.     }
  245.     for (unsigned short int i = 0; i <= n; i++)
  246.     {
  247.         massT[i] = yt(massX[i]);
  248.     }
  249.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  250.     cout << "Equation: y` = e^x + e^(-x); y(0) = 0; yt = e^x - e^(-x)" << endl;
  251.     cout << "Results: " << endl;
  252.     for (unsigned short int i = 0; i <= n; i++)
  253.     {
  254.         cout.precision(9);
  255.         cout << "X[" << i << "] = " << massX[i] << "\tY[" << i << "] = " << massY[i] << "\t"; printf("error = % 2.9f\n", fabs(massY[i] - massT[i]));
  256.     }
  257.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  258.     delete[] massX;
  259.     delete[] massY;
  260. }
  261. void milnMet()
  262. {
  263.     double a = 0.0, b = 1.0;
  264.     double n = (b - a) / h;
  265.     double* massX = new double[n + 1];
  266.     double* massY = new double[n + 1];
  267.     double* massT = new double[n + 1];
  268.     cout << "a (left border) = " << a << "\nb (right border) = " << b << endl;
  269.     cout << "h (step) = " << h << endl;
  270.     massX[0] = a;
  271.     for (unsigned short int i = 1; i <= n; i++)
  272.         massX[i] = a + i * h;
  273.  
  274.     massY[0] = 0.0; massY[1] = 0.200333507; massY[2] = 0.402672019; massY[3] = 0.609040608;
  275.     double temp;
  276.     for (unsigned short int i = 4; i <= n; i++)
  277.     {
  278.         temp = massY[i - 4] + (4 * h * (2 * derivative(massX[i - 3], massY[i - 3]) - derivative(massX[i - 2], massY[i - 2])
  279.             + 2 * derivative(massX[i - 1], massY[i - 1])) / 3.0);
  280.         massY[i] = massY[i - 2] + (h * (derivative(massX[i - 2], massY[i - 2]) + 4 * derivative(massX[i - 1], massY[i - 1]) +
  281.             derivative(massX[i], temp)) / 3.0);
  282.     }
  283.     for (unsigned short int i = 0; i <= n; i++)
  284.     {
  285.         massT[i] = yt(massX[i]);
  286.     }
  287.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  288.     cout << "Equation: y` = e^x + e^(-x); y(0) = 0; yt = e^x - e^(-x)" << endl;
  289.     cout << "Results: " << endl;
  290.     for (unsigned short int i = 0; i <= n; i++)
  291.     {
  292.         cout.precision(9);
  293.         cout << "X[" << i << "] = " << massX[i] << "\tY[" << i << "] = " << massY[i] << "\t"; printf("error = % 2.9f\n", fabs(massY[i] - massT[i]));
  294.     }
  295.     cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
  296.  
  297.     delete[] massX;
  298.     delete[] massY;
  299.     delete[] massT;
  300. }
  301. double rungeK(double h, double x, double y, int number)
  302. {
  303.     double result1, result2, result3, result4;
  304.     switch (number)
  305.     {
  306.     case 1:
  307.     {
  308.         result1 = h * derivative(x, y);
  309.         return result1;
  310.         break;
  311.     }
  312.     case 2:
  313.     {
  314.         result1 = h * derivative(x, y);
  315.         result2 = h * derivative(x + (h / 2.0), y + (result1 / 2.0));
  316.         return result2;
  317.         break;
  318.     }
  319.     case 3:
  320.     {
  321.         result1 = h * derivative(x, y);
  322.         result2 = h * derivative(x + (h / 2.0), y + (result1 / 2.0));
  323.         result3 = h * derivative(x + (h / 2.0), y + (result2 / 2.0));
  324.         return result3;
  325.         break;
  326.     }
  327.     case 4:
  328.     {
  329.         result1 = h * derivative(x, y);
  330.         result2 = h * derivative(x + (h / 2.0), y + (result1 / 2.0));
  331.         result3 = h * derivative(x + (h / 2.0), y + (result2 / 2.0));
  332.         result4 = h * derivative(x + h, y + result3);
  333.         return result4;
  334.         break;
  335.     }
  336.     default:cout << "ERROR_K" << endl;
  337.     }
  338. }
RAW Paste Data