API tools faq
paste
Khristina

Тест 4 "ур-е Гельмгольца" (внут. сфера) (новая формула)

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