Vladislav_Bezruk

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