Advertisement
Neveles

© 2021 Neveles. All rights reserved.

May 1st, 2021
948
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 3.42 KB | None | 0 0
  1. #include <iostream>
  2. #include <math.h>
  3. #include <iomanip>
  4.  
  5. #include "cmath.h"
  6. #include "quanc8.c"
  7. #include "zeroin.c"
  8. #include "rkf45.c"
  9.  
  10. double const u = 1 / 82.3;
  11.  
  12. double const A = 1.2;
  13. double const B = 0;
  14. double const C = 0;
  15. double D = 0;
  16.  
  17. double const T = 12;
  18. double const H = 0.4; // лучше брать меньше
  19.  
  20. double const ERR_FOR_QUANC = 1.0e-10;
  21. double const TOL_FOR_ZEROIN = 1.0e-10;
  22. double const ERR_FOR_RKF = 1.0e-4; // Больше брать нельзя
  23. double const MAXNFE = 30000;
  24.  
  25. double funForQuanc8(double z)
  26. {
  27.   return 1 / (exp(1.8 * z * z) + D);
  28. }
  29.  
  30. double quanc(double x)
  31. {
  32.   double a = 0.0;
  33.   double b = 1.0;
  34.   double abserr = ERR_FOR_QUANC;
  35.   double relerr = ERR_FOR_QUANC;
  36.   double result = 0.0;
  37.   double errest;
  38.   int nofun;
  39.   double posn;
  40.   int flag;
  41.  
  42.   double funForQuanc8(double z);
  43.  
  44.   D = x;
  45.   quanc8(funForQuanc8, a, b, abserr, relerr, &result, &errest, &nofun, &posn, &flag);
  46.  
  47.   return result - 0.3644632 * x;
  48. }
  49.  
  50. double calculateDWithZeroin()
  51. {
  52.   double d1 = 0.0;
  53.   double d2 = 3.0;
  54.   double tol = TOL_FOR_ZEROIN;
  55.   int flag;
  56.  
  57.   double quanc(double D);
  58.  
  59.   double resultD = zeroin(d1, d2, quanc, tol, &flag);
  60.  
  61.   std::cout.precision(8);
  62.   std::cout << "\nTOL = " << TOL_FOR_ZEROIN << '\n';
  63.   std::cout << "D   = " << std::setw(9) << D;
  64.   std::cout << "\nERR for rkf = " << ERR_FOR_RKF;
  65.   std::cout << "\nMAXNFE      = " << MAXNFE;
  66.   std::cout << "\n\n";
  67.  
  68.   return resultD;
  69. }
  70.  
  71. int funForRKF(int n, double t, double* x, double* dxdt)
  72. {
  73.   dxdt[0] = x[1]; // m1' = m2
  74.   dxdt[2] = x[3]; // m3' = m4
  75.  
  76.   dxdt[1] = -(1 - u) * x[0] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
  77.     (x[0] - 1) / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));
  78.  
  79.   dxdt[3] = -(1 - u) * x[2] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
  80.     x[2] / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));
  81.  
  82.   return 0;
  83. }
  84.  
  85. void RKF45(double* x0array, double* x1array, double* x2array, double* x3array)
  86. {
  87.   int n = 4;
  88.  
  89.   double x[4];
  90.   double xp[4];
  91.  
  92.   double t = 0.0;
  93.   double tout = 0.0;
  94.  
  95.   double relerr = ERR_FOR_RKF;
  96.   double abserr = ERR_FOR_RKF;
  97.  
  98.   double h = 0.0;
  99.  
  100.   int nfe = 0;
  101.   int maxnfe = 5000;
  102.  
  103.   int flag = 1;
  104.  
  105.   int step = 0;
  106.  
  107.   double D = calculateDWithZeroin();
  108.  
  109.   int  fail = 0;
  110.   rkfinit(n, &fail);
  111.  
  112.   if (fail != 0)
  113.   {
  114.     return;
  115.   }
  116.  
  117.   x[0] = A;
  118.   x[1] = B;
  119.   x[2] = C;
  120.   x[3] = D;
  121.  
  122.   std::cout << "   t        x          x'         y          y'\n";
  123.   std::cout << "------------------------------------------------------\n";
  124.   for (step = 1; step <= T / H; ++step)
  125.   {
  126.     tout = H * step;
  127.     t = tout - H;
  128.  
  129.     int funForRKF(int n, double t, double* x, double* dxdt);
  130.  
  131.     rkf45(funForRKF, n, x, xp, &t, tout, &relerr, abserr, &h, &nfe, maxnfe, &flag);
  132.  
  133.     printf("%5.1f %10.6f %10.6f %10.6f %10.6f \n", t, x[0], x[1], x[2], x[3]);
  134.  
  135.     if (flag != 2)
  136.     {
  137.       std::cout << "ERROR IN FLAG (flag = " << flag << ")\n";
  138.  
  139.       break;
  140.     }
  141.   }
  142.  
  143.   std::cout << "------------------------------------------------------\n";
  144.  
  145.  
  146.   rkfend();
  147.  
  148.   std::cout << "nfe  : " << nfe << '\n';
  149.   std::cout << "h    : " << h << '\n';
  150.   std::cout << "flag : " << flag << '\n';
  151.  
  152.   std::cout << "------------------------------------------------------\n\n";
  153. }
  154.  
  155. int main()
  156. {
  157.   double rkf45x0Array[30] = { 0 };
  158.   double rkf45x1Array[30] = { 0 };
  159.   double rkf45x2Array[30] = { 0 };
  160.   double rkf45x3Array[30] = { 0 };
  161.  
  162.   RKF45(rkf45x0Array, rkf45x1Array, rkf45x2Array, rkf45x3Array);
  163.  
  164.   return 0;
  165. }
  166.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement