• API
• FAQ
• Tools
• Archive
SHARE
TWEET

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

Khristina May 19th, 2019 81 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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. }
RAW Paste Data
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.

Top