Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <math.h>
- #include <iomanip>
- #include "cmath.h"
- #include "quanc8.c"
- #include "zeroin.c"
- #include "rkf45.c"
- double const u = 1 / 82.3;
- double const A = 1.2;
- double const B = 0;
- double const C = 0;
- double D = 0;
- double const T = 12;
- double const H = 0.4; // лучше брать меньше
- double const ERR_FOR_QUANC = 1.0e-10;
- double const TOL_FOR_ZEROIN = 1.0e-10;
- double const ERR_FOR_RKF = 1.0e-4; // Больше брать нельзя
- double const MAXNFE = 30000;
- double funForQuanc8(double z)
- {
- return 1 / (exp(1.8 * z * z) + D);
- }
- double quanc(double x)
- {
- double a = 0.0;
- double b = 1.0;
- double abserr = ERR_FOR_QUANC;
- double relerr = ERR_FOR_QUANC;
- double result = 0.0;
- double errest;
- int nofun;
- double posn;
- int flag;
- double funForQuanc8(double z);
- D = x;
- quanc8(funForQuanc8, a, b, abserr, relerr, &result, &errest, &nofun, &posn, &flag);
- return result - 0.3644632 * x;
- }
- double calculateDWithZeroin()
- {
- double d1 = 0.0;
- double d2 = 3.0;
- double tol = TOL_FOR_ZEROIN;
- int flag;
- double quanc(double D);
- double resultD = zeroin(d1, d2, quanc, tol, &flag);
- std::cout.precision(8);
- std::cout << "\nTOL = " << TOL_FOR_ZEROIN << '\n';
- std::cout << "D = " << std::setw(9) << D;
- std::cout << "\nERR for rkf = " << ERR_FOR_RKF;
- std::cout << "\nMAXNFE = " << MAXNFE;
- std::cout << "\n\n";
- return resultD;
- }
- int funForRKF(int n, double t, double* x, double* dxdt)
- {
- dxdt[0] = x[1]; // m1' = m2
- dxdt[2] = x[3]; // m3' = m4
- dxdt[1] = -(1 - u) * x[0] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
- (x[0] - 1) / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));
- dxdt[3] = -(1 - u) * x[2] / ((sqrt(x[0] * x[0] + x[2] * x[2])) * (x[0] * x[0] + x[2] * x[2])) -
- x[2] / (sqrt(pow(x[0] - 1, 2) + x[2] * x[2]) * (pow(x[0] - 1, 2) + x[2] * x[2]));
- return 0;
- }
- void RKF45(double* x0array, double* x1array, double* x2array, double* x3array)
- {
- int n = 4;
- double x[4];
- double xp[4];
- double t = 0.0;
- double tout = 0.0;
- double relerr = ERR_FOR_RKF;
- double abserr = ERR_FOR_RKF;
- double h = 0.0;
- int nfe = 0;
- int maxnfe = 5000;
- int flag = 1;
- int step = 0;
- double D = calculateDWithZeroin();
- int fail = 0;
- rkfinit(n, &fail);
- if (fail != 0)
- {
- return;
- }
- x[0] = A;
- x[1] = B;
- x[2] = C;
- x[3] = D;
- std::cout << " t x x' y y'\n";
- std::cout << "------------------------------------------------------\n";
- for (step = 1; step <= T / H; ++step)
- {
- tout = H * step;
- t = tout - H;
- int funForRKF(int n, double t, double* x, double* dxdt);
- rkf45(funForRKF, n, x, xp, &t, tout, &relerr, abserr, &h, &nfe, maxnfe, &flag);
- printf("%5.1f %10.6f %10.6f %10.6f %10.6f \n", t, x[0], x[1], x[2], x[3]);
- if (flag != 2)
- {
- std::cout << "ERROR IN FLAG (flag = " << flag << ")\n";
- break;
- }
- }
- std::cout << "------------------------------------------------------\n";
- rkfend();
- std::cout << "nfe : " << nfe << '\n';
- std::cout << "h : " << h << '\n';
- std::cout << "flag : " << flag << '\n';
- std::cout << "------------------------------------------------------\n\n";
- }
- int main()
- {
- double rkf45x0Array[30] = { 0 };
- double rkf45x1Array[30] = { 0 };
- double rkf45x2Array[30] = { 0 };
- double rkf45x3Array[30] = { 0 };
- RKF45(rkf45x0Array, rkf45x1Array, rkf45x2Array, rkf45x3Array);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement