• Sign Up
• Login
• API
• FAQ
• Tools
• Archive
Need a unique gift idea?
A Pastebin account makes a great Christmas gift
SHARE
TWEET

# Untitled

a guest Mar 14th, 2018 71 Never
ENDING IN00days00hours00mins00secs

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.

Top