API tools faq
paste
Khristina

Тест 1 "ур-е Лапласа" (внеш. сфера) (новая формула)

Khristina
Apr 12th, 2019
224
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 10.32 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <string>
  3. #include <cstring>
  4. #include <math.h>
  5. #include <cstdlib>
  6. #include <ctime>
  7. #include <cmath>
  8. #include <algorithm>
  9. using namespace std;
  10.  
  11. double epsilonJ(double H, double Q, double betta)
  12. {
  13.     double epsilon;
  14.     epsilon = (H / 2.0) + (Q * 1.0 / (betta * betta));
  15.     return epsilon;
  16. }
  17.  
  18. double betta0J(double delta, double H, double P, double betta)
  19. {
  20.     double betta0;
  21.     betta0 = ((delta * H) + (2 * P)) / (2.0 * betta * betta);
  22.     //cout << '\a';
  23.     return betta0;
  24. }
  25.  
  26. double hi0J(double H, double Q, double r, double betta)
  27. {
  28.     double hi0;
  29.     hi0 = (H * H / 4.0) + (((H * Q) + (r * r)) / (betta * betta));
  30.     return hi0;
  31. }
  32.  
  33. double I3J(double epsilon1, double deltaplus, double deltaminus, double zplus, double zminus)
  34. {
  35.     double I3 = 0;
  36.     if (fabs(deltaplus) > pow(10, -10))
  37.     {
  38.         if (((epsilon1 * epsilon1) - (deltaplus * deltaminus)) > pow(10, -10)) //+
  39.         {
  40.             I3 = (log(fabs(((deltaplus * zplus) + epsilon1 - sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus))) / (((deltaplus * zplus) + epsilon1 + sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus)))))) / (2 * sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus)))) - (log(fabs(((deltaplus * zminus) + epsilon1 - sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus))) / (((deltaplus * zminus) + epsilon1 + sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus)))))) / (2 * sqrt((epsilon1 * epsilon1) - (deltaplus * deltaminus))));
  41.         }
  42.  
  43.         else if (((epsilon1 * epsilon1) - (deltaplus * deltaminus)) < -pow(10, -10))
  44.         {
  45.             I3 = (atan(((deltaplus * zplus) + epsilon1) / (sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1)))) / sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1))) - (atan(((deltaplus * zminus) + epsilon1) / (sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1)))) / sqrt((deltaplus * deltaminus) - (epsilon1 * epsilon1)));
  46.         }
  47.  
  48.         else if (fabs((epsilon1 * epsilon1) - (deltaplus * deltaminus)) <= pow(10, -10))
  49.         {
  50.             I3 = ((-1.0) / ((deltaplus * zplus) + epsilon1)) + (1.0 / (deltaplus * zminus + epsilon1));
  51.         }
  52.     }
  53.     else if (fabs(deltaplus) <= pow(10, -10))
  54.     {
  55.         if (fabs(epsilon1) > pow(10, -10))
  56.         {
  57.             I3 = (log(fabs((2 * epsilon1 * zplus) + deltaminus)) / (2.0 * epsilon1)) - (log(fabs((2 * epsilon1 * zminus) + deltaminus)) / (2 * epsilon1));
  58.         }
  59.         else
  60.         {
  61.             I3 = (zplus / deltaminus) - (zminus / deltaminus);
  62.         }
  63.     }
  64.  
  65.     return I3;
  66. }
  67.  
  68. double I2(double deltaplus, double zplus, double zminus, double epsilon1, double deltaminus)
  69. {
  70.     double I;
  71.     I = I3J(epsilon1, deltaplus, deltaminus, zplus, zminus);
  72.     double I2;
  73.     I2 = (((zplus - (epsilon1 * log(fabs((deltaplus * zplus * zplus) + (2 * epsilon1 * zplus) + deltaminus)) / deltaplus)) - (zminus - (epsilon1 * log(fabs((deltaplus * zminus * zminus) + (2 * epsilon1 * zminus) + deltaminus)) / deltaplus))) / deltaplus) + ((2 * epsilon1 * epsilon1) - (deltaplus * deltaminus)) * I / (deltaplus * deltaplus);
  74.     return I2;
  75. }
  76.  
  77. double I4(double deltaminus, double zplus, double zminus, double epsilon1, double deltaplus)
  78. {
  79.     double I4;
  80.     double I;
  81.     I = I3J(epsilon1, deltaplus, deltaminus, zplus, zminus);
  82.     if (fabs(deltaplus) > pow(10, -10))
  83.     {
  84.         I4 = ((((-1.0) / zplus) - (epsilon1 * log(zplus * zplus / fabs((deltaplus * zplus * zplus) + (2 * epsilon1 * zplus) + deltaminus)) / deltaminus) - ((-1.0) / zminus) + (epsilon1 * log(zminus * zminus / fabs((deltaplus * zminus * zminus) + (2 * epsilon1 * zminus) + deltaminus)) / deltaminus)) / deltaminus) + (((2 * epsilon1 * epsilon1) - (deltaplus * deltaminus)) * I / (deltaminus * deltaminus));
  85.     }
  86.     else
  87.     {
  88.         if (fabs(epsilon1) > pow(10, -10))
  89.         {
  90.             I4 = (((-1) / (deltaminus * zplus)) + ((2 * epsilon1 * log(fabs(((2 * epsilon1 * zplus) + deltaminus) / zplus))) / (deltaminus * deltaminus))) - (((-1) / (deltaminus * zminus)) + (2 * epsilon1 * log(fabs((2 * epsilon1 * zminus + deltaminus) / zminus))) / (deltaminus * deltaminus));
  91.         }
  92.  
  93.         else
  94.         {
  95.             I4 = ((-1.0) / (deltaminus * zplus)) + (1.0 / (deltaminus * zminus));
  96.         }
  97.     }
  98.  
  99.     return I4;
  100. }
  101.  
  102. double I1(double deltaplus, double hi1, double deltaminus, double zplus, double zminus, double epsilon1)
  103. {
  104.     double I1;
  105.     double II, III, IIII;
  106.     II = I2(deltaplus, zplus, zminus, epsilon1, deltaminus);
  107.     III = I3J(epsilon1, deltaplus, deltaminus, zplus, zminus);
  108.     IIII = I4(deltaminus, zplus, zminus, epsilon1, deltaplus);
  109.     I1 = (deltaplus * II / 2.0) - (hi1 * III) + (deltaminus * IIII / 2.0);
  110.     return I1;
  111. }
  112.  
  113. double I5(double U1, double U2, double k1, double k2)
  114. {
  115.     double I5;
  116.     if (fabs(k1) > pow(10, -10))
  117.     {
  118.         I5 = (-U2 + ((U2 + (k2 / k1)) * log(fabs(k1 * U2 + k2)))) - (-U1 + ((U1 + (k2 / k1)) * log(fabs(k1 * U1 + k2))));
  119.     }
  120.     else
  121.     {
  122.         I5 = (U2 - U1) * log(fabs(k2));
  123.     }
  124.     return I5;
  125. }
  126.  
  127. int main()
  128. {
  129.     int M = 25, N = 25;
  130.     long double V0 = 0, S0 = 0;
  131.     long double p = 0;
  132.     const double pi = 3.14159265358979;
  133.     double x1, x2, x3, x, y1, y2, y3, y11, y21, y31, h = 2.0 * pi / N, H = pi / (1.0 * M), etta, R;
  134.     double A = N * h, B = M * H;
  135.     double tetta, betta, y1prv, y2prv, y3prv, alfa, y1pru, y2pru, y3pru, delta;
  136.     double P, Q, r1, r2, r3, r;
  137.     double delta0, alfa0, epsilon, betta0, hi0, hi1;
  138.     double deltaplus, deltaminus, delta1, epsilon1;
  139.     double IH = 0, zminus, zplus, IminusH = 0;
  140.     double dzettaminus, dzettaplus;
  141.     setlocale(LC_ALL, "Rus");
  142.     for (int i = 1; i < 6; i++)
  143.     {
  144.         R = 1 + pow(10, (-1) * i);
  145.         cout << "R=" << R << endl;
  146.         p = 0;
  147.         for (int q = 0; q <= 2; q++)
  148.         {
  149.             double uq = pi * q / (1.0 * N);
  150.  
  151.             for (int l = 0; l <= 2 * M; l++)
  152.             {
  153.                 double vl = pi * l / (2.0 * M);
  154.                 y11 = sin(vl) * cos(uq); // !
  155.                 y21 = sin(vl) * sin(uq); // !
  156.                 y31 = cos(vl); // !
  157.                 x1 = y11 * R;
  158.                 x2 = y21 * R;
  159.                 x3 = y31 * R;
  160.                 x = sqrt((x1 * x1) + (x2 * x2) + (x3 * x3));
  161.  
  162.                 for (int m = 0; m < M; m++)
  163.                 {
  164.                     double vm = (m + 0.5) * H;
  165.                     for (int n = 0; n < N; n++)
  166.                     {
  167.                         double un = (n + 0.5) * h;
  168.                         y1 = sin(vm) * cos(un);
  169.                         y2 = sin(vm) * sin(un);
  170.                         y3 = cos(vm);
  171.                         y1prv = cos(un) * cos(vm);
  172.                         y2prv = sin(un) * cos(vm);
  173.                         y3prv = (-1) * sin(vm);
  174.                         y1pru = (-1) * sin(un) * sin(vm);
  175.                         y2pru = sin(vm) * cos(un);
  176.                         y3pru = 0;
  177.                         r1 = y1 - x1;
  178.                         r2 = y2 - x2;
  179.                         r3 = y3 - x3;
  180.                         r = sqrt((r1 * r1) + (r2 * r2) + (r3 * r3));
  181.                         betta = sqrt((y1prv * y1prv) + (y2prv * y2prv) + (y3prv * y3prv));
  182.                         etta = sin(vm);
  183.                         alfa = sqrt((y1pru * y1pru) + (y2pru * y2pru) + (y3pru * y3pru));
  184.                         delta = (y1pru * y1prv) + (y2pru * y2prv) + (y3pru * y3prv);
  185.                         P = (r1 * y1pru) + (r2 * y2pru) + (r3 * y3pru);
  186.                         Q = (r1 * y1prv) + (r2 * y2prv) + (r3 * y3prv);
  187.                         delta0 = delta * 1.0 / (betta * betta);
  188.                         alfa0 = fabs(alfa / betta); //!
  189.                         epsilon = epsilonJ(H, Q, betta);
  190.                         betta0 = betta0J(delta, H, P, betta);
  191.                         hi0 = hi0J(H, Q, r, betta);
  192.                         hi1 = sqrt(hi0 - (betta0 * betta0 / (alfa0 * alfa0)));
  193.                         epsilon1 = epsilon - (delta0 * betta0 / (alfa0 * alfa0));
  194.                         delta1 = delta0 * hi1 / alfa0;
  195.                         deltaplus = hi1 + delta1;
  196.                         deltaminus = hi1 - delta1;
  197.                         dzettaplus = ((alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
  198.                         dzettaminus = (((-1) * alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
  199.                         zplus = dzettaplus + sqrt((dzettaplus * dzettaplus) + 1);
  200.                         zminus = dzettaminus + sqrt((dzettaminus * dzettaminus) + 1);
  201.                         if (hi1 > pow(10, -10))
  202.                         {
  203.                             IH = hi1 * ((((dzettaplus * log(fabs(epsilon1 + (delta1 * dzettaplus) + (hi1 * sqrt((dzettaplus * dzettaplus) + 1))))) - (dzettaminus * log(fabs(epsilon1 + (delta1 * dzettaminus) + (hi1 * sqrt((dzettaminus * dzettaminus) + 1))))))) - I1(deltaplus, hi1, deltaminus, zplus, zminus, epsilon1)) / alfa0;
  204.                         }
  205.                         else if (fabs(hi1) <= pow(10, -10))
  206.                         {
  207.                             if (-h / 2 < (-betta0 / (alfa0 * alfa0)) && (-betta0 / (alfa0 * alfa0)) < h / 2)
  208.                             {
  209.                                 IH = I5(-h / 2, (-betta0 / (alfa0 * alfa0)), delta0 - alfa0, epsilon - (betta0 / alfa0)) + I5((-betta0 / (alfa0 * alfa0)), h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
  210.                             }
  211.                             else if ((-betta0 / (alfa0 * alfa0)) >= h / 2)
  212.                             {
  213.                                 IH = I5(-h / 2, h / 2, delta0 - alfa0, epsilon - (betta0 / alfa0));
  214.                             }
  215.                             else if ((-betta0 / (alfa0 * alfa0)) <= -h / 2)
  216.                             {
  217.                                 IH = I5(-h / 2, h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
  218.                             }
  219.                         }
  220.                         epsilon = epsilonJ(-H, Q, betta);
  221.                         betta0 = betta0J(delta, -H, P, betta);
  222.                         hi0 = hi0J(-H, Q, r, betta);
  223.                         hi1 = sqrt(hi0 - (betta0 * betta0 / (alfa0 * alfa0)));
  224.                         epsilon1 = epsilon - (delta0 * betta0 / (alfa0 * alfa0));
  225.                         delta1 = delta0 * hi1 / alfa0;
  226.                         deltaplus = hi1 + delta1;
  227.                         deltaminus = hi1 - delta1;
  228.                         dzettaplus = ((alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
  229.                         dzettaminus = (((-1) * alfa0 * h / 2) + (betta0 / alfa0)) / hi1;
  230.                         zplus = dzettaplus + sqrt(dzettaplus * dzettaplus + 1);
  231.                         zminus = dzettaminus + sqrt(dzettaminus * dzettaminus + 1);
  232.                         if (hi1 > pow(10, -10))
  233.                         {
  234.                             IminusH = hi1 * ((((dzettaplus * log(fabs(epsilon1 + (delta1 * dzettaplus) + (hi1 * sqrt((dzettaplus * dzettaplus) + 1))))) - (dzettaminus * log(fabs(epsilon1 + (delta1 * dzettaminus) + (hi1 * sqrt((dzettaminus * dzettaminus) + 1))))))) - I1(deltaplus, hi1, deltaminus, zplus, zminus, epsilon1)) / alfa0;
  235.                         }
  236.                         else if (fabs(hi1) <= pow(10, -10))
  237.                         {
  238.                             if (-h / 2 < (-betta0 / (alfa0 * alfa0)) && (-betta0 / (alfa0 * alfa0)) < h / 2)
  239.                             {
  240.                                 IminusH = I5(-h / 2, (-betta0 / (alfa0 * alfa0)), delta0 - alfa0, epsilon - (betta0 / alfa0)) + I5((-betta0 / (alfa0 * alfa0)), h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
  241.                             }
  242.                             else if ((-betta0 / (alfa0 * alfa0)) >= h / 2)
  243.                             {
  244.                                 IminusH = I5(-h / 2, h / 2, delta0 - alfa0, epsilon - (betta0 / alfa0));
  245.                             }
  246.                             else if ((-betta0 / (alfa0 * alfa0)) <= -h / 2)
  247.                             {
  248.                                 IminusH = I5(-h / 2, h / 2, delta0 + alfa0, epsilon + (betta0 / alfa0));
  249.                             }
  250.                         }
  251.                         tetta = 1.0 * (IH - IminusH) / betta;
  252.                         S0 = S0 + (etta * tetta);
  253.                     }
  254.                 }
  255.                 if (fabs(x) < 1)
  256.                 {
  257.                     V0 = 4 * pi;
  258.                 }
  259.                 else if (fabs(x) > 1)
  260.                 {
  261.                     V0 = 4.0 * pi / fabs(x);
  262.                 }
  263.                 //cout << S0 << endl;
  264.                 //cout << V0 << endl;
  265.                 long double maximum = max(p, fabs(S0 - V0) * 1.0 / V0);
  266.                 S0 = 0;
  267.                 p = maximum;
  268.             }
  269.         }
  270.         cout << "Относительная погрешность=" << p << endl;
  271.         cout << "______________________" << endl;
  272.     }
  273.  
  274.     return 0;
  275. }
Add Comment
Please, Sign In to add comment
Public Pastes
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.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!