# Runge–Kutta method

Nov 29th, 2021 (edited)
810
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #include <iostream>
2. #include <cmath>
3.
4. using namespace std;
5.
6. struct Point {
7.     double x;
8.     double y;
9.
10.     Point() {x = y = 0;}
11.     Point(double _x, double _y) : x(_x), y(_y) {};
12. };
13.
14. double df(Point x) {return 0.2 * x.x + pow(x.y, 2);}
15.
16. void print(Point* x, int n) {
17.     cout << endl << "n = " << n << endl << endl;
18.     cout << "===== points =====" << endl;
19.     for (int i = 0; i < n + 1; i++) {
20.         cout << "[" << i << "]: " << "x = " << x[i].x << " y = " << x[i].y << endl;
21.         cout << endl;
22.     }
23. }
24.
25. Point* rkSolve(Point s, double a, double b, double df(Point x), int n) {
26.     double k1, k2, k3, k4;
27.     double h = (b - a) / n;
28.     Point* res = new Point[n + 1];
29.
30.     for (int i = 0; i < n + 1; i++) {
31.         res[i] = s;
32.
33.         k1 = df(s);
34.         k2 = df(Point(s.x + 0.5 * h, s.y + 0.5 * h * k1));
35.         k3 = df(Point(s.x + 0.5 * h, s.y + 0.5 * h * k2));
36.         k4 = df(Point(s.x + h, s.y + h * k3));
37.
38.         s.y = s.y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
39.         s.x += h;
40.     }
41.     return res;
42. }
43.
44. double checkEps(Point* s1, Point* s2, int n) {
45.     double eps = -1, curEps;
46.
47.     for (int i = 1; i < n + 1; i++) {
48.         curEps = fabs(s1[i].y - s2[i * 2].y) / 15;
49.         eps = max(eps, curEps);
50.     }
51.     return eps;
52. }
53.
54. void solve(Point s, double a, double b, double df(Point x), double eps) {
55.     int n = 10;
56.     Point *s1, *s2;
57.
58.     do {
59.         n *= 2;
60.
61.         s1 = rkSolve(s, a, b, df, n);
62.         s2 = rkSolve(s, a, b, df, n * 2);
63.
64.     } while (checkEps(s1, s2, n) >= eps);
65.
66.     print(s1, n);
67. }
68.
69. int main() {
70.     Point s;
71.     double a, b, eps;
72.
73.     cout << "Enter start point: ";
74.     cin >> s.x >> s.y;
75.
76.     cout << "Enter a & b: ";
77.     cin >> a >> b;
78.
79.     cout << "Enter eps: ";
80.     cin >> eps;
81.
82.     solve(s, a, b, df, eps);
83.
84.     return 0;
85. }
RAW Paste Data