daily pastebin goal
19%
SHARE
TWEET

Untitled

a guest Mar 14th, 2018 71 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <algorithm>
  3. #include <cmath>
  4. #include <functional>
  5.  
  6. using namespace std;
  7.  
  8. function <double(double)> y = [](double x) {return 3.0 / 2.0 * exp(-x) + x / 2.0 - 1.0 / 2.0;};
  9. function <double(double)> dy = [](double x) {return 1.0 / 2.0 - 3.0 / 2.0 * exp(-x);};
  10. function <double(double)> ddy = [](double x) {return 3.0 / 2.0 * exp(-x);};
  11. function <double(double)> dddy = [](double x) {return -3.0 / 2.0 * exp(-x);};
  12. function <double(double)> ddddy = ddy;
  13. function <double(double)> dddddy = dddy;
  14. function <double(double)> ddddddy = ddddy;
  15. function <double(double, double)> f = [](double x, double y) {return -y + x / 2.0;};
  16.  
  17. struct Point {
  18.     double node;
  19.     double value;
  20. };
  21.  
  22. struct DistanceFromPointToX
  23. {
  24.     double distance;
  25.     int pointIndex;
  26. };
  27.  
  28. bool cmp(DistanceFromPointToX a, DistanceFromPointToX b) {
  29.     return a.distance < b.distance;
  30. }
  31.  
  32. int factorial (int n) {
  33.     int result = 1;
  34.     for (int i = 1; i <= n; ++i)
  35.         result *= i;
  36.     return result;
  37. }
  38.  
  39. vector <Point> orderPoints(const double x,
  40.                            vector <Point> points) {
  41.     vector <DistanceFromPointToX> distances;
  42.     for (int i = 0; i < points.size(); ++i)
  43.         distances.push_back({fabs(points[i].node - x), i});
  44.     sort(distances.begin(), distances.end(), cmp);
  45.     vector <Point> sortedPoints;
  46.     for (int i = 0; i < distances.size(); ++i)
  47.         sortedPoints.push_back(points[distances[i].pointIndex]);
  48.     return sortedPoints;
  49. }
  50.  
  51. void getFiniteDifferenceTable(vector <Point> points,
  52.                               const int n,
  53.                               vector <double> *coefficients,
  54.                               double h) {
  55.     for (int i = 0; i <= n; ++i)
  56.         coefficients[0].push_back(points[i].value * h);
  57.     for (int i = 1; i <= n; ++i)
  58.         for (int j = 0; j <= n - i; ++j)
  59.             coefficients[i].push_back(coefficients[i - 1][j + 1] - coefficients[i - 1][j]);
  60. }
  61.  
  62. double AdamsRule(const double x,
  63.                  vector <Point> points,
  64.                  const int n,
  65.                  const double h,
  66.                  vector <double> ys) {
  67.     vector <double> coefficients[n + 1];
  68.     getFiniteDifferenceTable(points, n, coefficients, h);
  69.     vector <double> etas;
  70.     for (int i = 0; i <= 4; ++i)
  71.         etas.push_back(coefficients[i][n - i]);
  72.     double result = ys[ys.size() - 1] + etas[0] + etas[1] * 1.0/2.0 + etas[2] * 5.0/12.0 + etas[3] * 3.0/8.0 + etas[4] * 251.0/720.0;
  73.     return result;
  74. }
  75.  
  76. void FirstMeth(const int n, const double x_0, const double h, const double y_0) {
  77.     vector <double> ys;
  78.     ys.clear();
  79.     ys.push_back(y_0);
  80.     for (int i = 1; i <= n; ++i) {
  81.         double x = x_0 + (i - 1) * h;
  82.         double currentY = y(x) + h / 2.0 * f(x, y(x));
  83.         ys.push_back(y(x) + h * f(x + h / 2.0, currentY));
  84.     }
  85.  
  86.     for (int i = 0; i < ys.size(); ++i) {
  87.         double x = x_0 + i * h;
  88.         cout << "y(" << x << ") = " << ys[i] << ", абсолютная погрешность " << fabs(y(x) - ys[i]) << endl;
  89.     }
  90.     cout << endl << endl << endl;
  91. }
  92.  
  93. void SecondMeth(const int n, const double x_0, const double h, const double y_0) {
  94.     vector <double> ys, cauchyY;
  95.     ys.clear(), cauchyY.clear();
  96.     ys.push_back(y_0), cauchyY.push_back(y_0);
  97.     for (int i = 1; i <= n; ++i) {
  98.         double x = x_0 + (i - 1) * h;
  99.         ys.push_back(ys[i - 1] + h * f(x - h, ys[i - 1]));
  100.     }
  101.     for (int i = 1; i <= n; ++i) {
  102.         double x = x_0 + (i - 1) * h;
  103.         cauchyY.push_back(cauchyY[i - 1] + h / 2.0 * (f(x - h, cauchyY[i - 1]) + f(x, ys[i])));
  104.     }
  105.  
  106.     for (int i = 0; i < cauchyY.size(); ++i) {
  107.         double x = x_0 + i * h;
  108.         cout << "y(" << x << ") = " << cauchyY[i] << ", абсолютная погрешность " << fabs(y(x) - cauchyY[i]) << endl;
  109.     }
  110.     cout << endl << endl << endl;
  111. }
  112.  
  113.  
  114. int main() {
  115.  
  116.     setlocale(LC_ALL, "Russian");
  117.     cout << "Численное решение Задачи Коши для обыкновенного дифференциального уравнения первого порядка" << endl;
  118.  
  119.     cout << "Дифференциальное уравнение: y'(x) = -y(x) + x / 2" << endl;
  120.     cout << "Начальное условие: y(0) = 1" << endl;
  121.     cout << "Шаг: h = 0.1" << endl;
  122.     cout << "Количество узлов: N = 10" << endl << endl;
  123.  
  124.     cout << "Решение дифференциального уравнения в общем виде: y(x) = 3 / 2 * exp(-x) + x / 2 - 1 / 2." << endl;
  125.  
  126.     double x_0 = 0.0;
  127.     double y_0 = 1.0;
  128.     double h = 0.1;
  129.     const int n = 10;
  130.  
  131.     cout << "Точное решение: " << endl;
  132.     vector <double> ys;
  133.     for (int i = -2; i <= n; ++i) {
  134.         if (i <= 2)
  135.             ys.push_back(y(x_0 + i * h));
  136.         cout << "y(" << x_0 + (double)i * h << ") = " << y(x_0 + (double)i * h) << endl;
  137.     }
  138.     cout << endl << endl << endl;
  139.  
  140.     cout << "Метод разложения в ряд Тейлора: " << endl << endl;
  141.  
  142.     for (int i = -2; i <= 2; ++i){
  143.         double x = x_0 + i * h;
  144.         double curResult = y(x_0) * pow(x - x_0, 0) / (double)factorial(0) +
  145.                 dy(x_0) * pow(x - x_0, 1) / (double)factorial(1) +
  146.                 ddy(x_0) * pow(x - x_0, 2) / (double)factorial(2) +
  147.                 dddy(x_0) * pow(x - x_0, 3) / (double)factorial(3) +
  148.                 ddddy(x_0) * pow(x - x_0, 4) / (double)factorial(4) +
  149.                 dddddy(x_0) * pow(x - x_0, 5) / (double)factorial(5) +
  150.                 ddddddy(x_0) * pow(x - x_0, 6) / (double)factorial(6);
  151.         cout << "y(" << x << ") = " << curResult << ", абсолютная погрешность " << fabs(y(x) - curResult) << endl;
  152.     }
  153.     cout << endl << endl << endl;
  154.  
  155.  
  156.     cout << "Экстраполяционный метод Адамса 4-го порядка: " << endl;
  157.     vector <Point> points;
  158.     for (int i = -2; i <= 2; ++i) {
  159.         double x = x_0 + i * h;
  160.         points.push_back({x, -y(x) + x / 2.0});
  161.     }
  162.     for (int i = 3; i <= n; ++i) {
  163.         double x = x_0 + i * h;
  164.         double result = AdamsRule(x, points, i + 1, h, ys);
  165.         ys.push_back(result);
  166.         points.push_back({x, -result + x / 2.0});
  167.         cout << "y(" << x << ") = " << result << ", абсолютная погрешность " << fabs(y(x) - result) << endl;
  168.     }
  169.     cout << endl << endl << endl;
  170.  
  171.  
  172.     cout << "Метод Рунге-Кутта 4-го порядка: " << endl;
  173.     ys.clear();
  174.     ys.push_back(y_0);
  175.     for (int i = 1; i <= n; ++i) {
  176.         double x = x_0 + (i - 1) * h;
  177.         double k1 = h * f(x, y(x));
  178.         double k2 = h * f(x + h / 2.0, y(x) + k1 / 2.0);
  179.         double k3 = h * f(x + h / 2.0, y(x) + k2 / 2.0);
  180.         double k4 = h * f(x + h, y(x) + k3);
  181.         ys.push_back(ys[ys.size() - 1] + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0);
  182.     }
  183.     for (int i = 0; i < ys.size(); ++i) {
  184.         double x = x_0 + i * h;
  185.         cout << "y(" << x << ") = " << ys[i] << ", абсолютная погрешность " << fabs(y(x) - ys[i]) << endl;
  186.     }
  187.     cout << endl << endl << endl;
  188.  
  189.  
  190.     cout << "Методы Эйлера: " << endl;
  191.     cout << "Нулевой метод: " << endl;
  192.     for (int i = 1; i <= n; ++i) {
  193.         double x = x_0 + (i - 1) * h;
  194.         double result = y(x) + h * f(x, y(x));
  195.         cout << "y(" << x + h << ") = " << result << ", абсолютная погрешность " << fabs(y(x) - result) << endl;
  196.     }
  197.     cout << endl << endl << endl;
  198.  
  199.     cout << "Первый метод: " << endl;
  200.     FirstMeth(n, x_0, h, y_0);
  201.  
  202.  
  203.     vector <double> ys1;
  204.  
  205.     ys1.clear();
  206.     cout << "Второй метод: " << endl;
  207.     SecondMeth(n, x_0, h, y_0);
  208.  
  209.     return 0;
  210. }
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