Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <string>
- #include <cstring>
- #include <math.h>
- #include <cstdlib>
- #include <ctime>
- #include <cmath>
- #include <complex>
- #include <algorithm>
- using namespace std;
- double epsilonJ(double H, double Q, double betta)
- {
- double epsilon;
- epsilon = (H / 2.0) + (Q * 1.0 / (betta * betta));
- return epsilon;
- }
- double betta0J(double delta, double H, double P, double betta)
- {
- double betta0;
- betta0 = ((delta * H) + (2 * P)) / (2.0 * betta * betta);
- //cout << '\a';
- return betta0;
- }
- double hi0J(double H, double Q, double r, double betta)
- {
- double hi0;
- hi0 = (H * H / 4.0) + (((H * Q) + (r * r)) / (betta * betta));
- return hi0;
- }
- double I3J(double epsilon1, double deltaplus, double deltaminus, double zplus, double zminus)
- {
- double I3 = 0;
- if (fabs(deltaplus) > pow(10, -10))
- {
- if (((epsilon1 * epsilon1) - (deltaplus * deltaminus)) > pow(10, -10)) //+
- {
- 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))));
- }
- else if (((epsilon1 * epsilon1) - (deltaplus * deltaminus)) < -pow(10, -10))
- {
- 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)));
- }
- else if (fabs((epsilon1 * epsilon1) - (deltaplus * deltaminus)) <= pow(10, -10))
- {
- I3 = (1.0 / ((deltaplus * zminus) + epsilon1)) - (1.0 / ((deltaplus * zplus) + epsilon1));
- }
- }
- else if (fabs(deltaplus) <= pow(10, -10))
- {
- if (fabs(epsilon1) > pow(10, -10))
- {
- I3 = (log(fabs((2 * epsilon1 * zplus) + deltaminus)) / (2.0 * epsilon1)) - (log(fabs((2.0 * epsilon1 * zminus) + deltaminus)) / (2.0 * epsilon1));
- }
- else
- {
- I3 = (zplus / deltaminus) - (zminus / deltaminus);
- }
- }
- return I3;
- }
- double I2(double deltaplus, double zplus, double zminus, double epsilon1, double deltaminus)
- {
- double I;
- I = I3J(epsilon1, deltaplus, deltaminus, zplus, zminus);
- double I2;
- 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);
- if (fabs(deltaplus) < pow(10, -10))
- {
- I2 = 0;
- }
- return I2;
- }
- double I4(double deltaminus, double zplus, double zminus, double epsilon1, double deltaplus)
- {
- double I4;
- double I;
- I = I3J(epsilon1, deltaplus, deltaminus, zplus, zminus);
- if (fabs(deltaplus) > pow(10, -10))
- {
- 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));
- }
- else
- {
- if (fabs(epsilon1) > pow(10, -10))
- {
- 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);
- }
- else
- {
- I4 = (1.0 / (deltaminus * zminus)) - (1.0 / (deltaminus * zplus));
- }
- }
- return I4;
- }
- double I1(double deltaplus, double hi1, double deltaminus, double zplus, double zminus, double epsilon1)
- {
- double I1;
- double II, III, IIII;
- II = I2(deltaplus, zplus, zminus, epsilon1, deltaminus);
- III = I3J(epsilon1, deltaplus, deltaminus, zplus, zminus);
- IIII = I4(deltaminus, zplus, zminus, epsilon1, deltaplus);
- I1 = (deltaplus * II / 2.0) - (hi1 * III) + (deltaminus * IIII / 2.0);
- return I1;
- }
- double I5(double U1, double U2, double k1, double k2)
- {
- double I5;
- if (fabs(k1) > pow(10, -10))
- {
- I5 = (-U2 + ((U2 + (1.0 * k2 / k1)) * log(fabs((k1 * U2) + k2)))) + U1 - ((U1 + (1.0 * k2 / k1)) * log(fabs((k1 * U1) + k2)));
- }
- else
- {
- I5 = (U2 - U1) * log(fabs(k2));
- }
- return I5;
- }
- int main()
- {
- int M = 25, N = 25;
- //long double S0 = 0;
- //long double p = 0;
- const double pi = 3.14159265358979;
- double x1, x2, x3, x, y1, y2, y3, y11, y21, y31, h = 2.0 * pi / N, H = pi / (1.0 * M), etta, R, xminusy;
- double A = N * h, B = M * H;
- double tetta, betta, y1prv, y2prv, y3prv, alfa, y1pru, y2pru, y3pru, delta;
- double P, Q, r1, r2, r3, r;
- double delta0, alfa0, epsilon, betta0, hi0, hi1;
- double deltaplus, deltaminus, delta1, epsilon1;
- double IH = 0, zminus, zplus, IminusH = 0;
- double dzettaminus, dzettaplus;
- complex< double> V0(0, 0);
- complex< double> S0(0, 0);
- complex< double > z(0, 1.0); // мнимая единица;
- complex <double> p(0, 0);
- setlocale(LC_ALL, "Rus");
- for (int i = 1; i < 6; i++)
- {
- R = 1 + pow(10, (-1) * i);
- cout << "R=" << R << endl;
- p = 0;
- for (int q = 0; q <= 2; q++)
- {
- double uq = pi * q / (1.0 * N);
- for (int l = 0; l <= 2 * M; l++)
- {
- double vl = pi * l / (2.0 * M);
- y11 = sin(vl) * cos(uq); // !
- y21 = sin(vl) * sin(uq); // !
- y31 = cos(vl); // !
- x1 = y11 * R;
- x2 = y21 * R;
- x3 = y31 * R;
- x = sqrt((x1 * x1) + (x2 * x2) + (x3 * x3));
- for (int m = 0; m < M; m++)
- {
- double vm = (m + 0.5) * H;
- for (int n = 0; n < N; n++)
- {
- double un = (n + 0.5) * h;
- y1 = sin(vm) * cos(un);
- y2 = sin(vm) * sin(un);
- y3 = cos(vm);
- xminusy = pow(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2), 0.5);
- y1prv = cos(un) * cos(vm);
- y2prv = sin(un) * cos(vm);
- y3prv = (-1) * sin(vm);
- y1pru = (-1) * sin(un) * sin(vm);
- y2pru = sin(vm) * cos(un);
- y3pru = 0.0;
- r1 = y1 - x1;
- r2 = y2 - x2;
- r3 = y3 - x3;
- r = sqrt((r1 * r1) + (r2 * r2) + (r3 * r3));
- betta = sqrt((y1prv * y1prv) + (y2prv * y2prv) + (y3prv * y3prv));
- etta = sin(vm);
- alfa = sqrt((y1pru * y1pru) + (y2pru * y2pru) + (y3pru * y3pru));
- delta = (y1pru * y1prv) + (y2pru * y2prv) + (y3pru * y3prv);
- P = (r1 * y1pru) + (r2 * y2pru) + (r3 * y3pru);
- Q = (r1 * y1prv) + (r2 * y2prv) + (r3 * y3prv);
- delta0 = delta * 1.0 / (betta * betta);
- alfa0 = fabs(alfa / betta); //!
- epsilon = epsilonJ(H, Q, betta);
- betta0 = betta0J(delta, H, P, betta);
- hi0 = hi0J(H, Q, r, betta);
- hi1 = sqrt(hi0 - (betta0 * betta0 / (alfa0 * alfa0)));
- epsilon1 = epsilon - (delta0 * betta0 / (alfa0 * alfa0));
- delta1 = delta0 * hi1 / alfa0;
- deltaplus = hi1 + delta1;
- deltaminus = hi1 - delta1;
- dzettaplus = ((alfa0 * h / 2.0) + (betta0 / alfa0)) / hi1;
- dzettaminus = (((-1) * alfa0 * h / 2.0) + (betta0 / alfa0)) / hi1;
- zplus = dzettaplus + sqrt((dzettaplus * dzettaplus) + 1);
- zminus = dzettaminus + sqrt((dzettaminus * dzettaminus) + 1);
- if (hi1 > pow(10, -10))
- {
- 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;
- }
- else if (fabs(hi1) <= pow(10, -10))
- {
- if (-h / 2.0 < (-betta0 / (alfa0 * alfa0)) && (-betta0 / (alfa0 * alfa0)) < h / 2.0)
- {
- 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));
- }
- else if ((-betta0 / (alfa0 * alfa0)) >= h / 2)
- {
- IH = I5(-h / 2.0, h / 2.0, delta0 - alfa0, epsilon - (betta0 / alfa0));
- }
- else if ((-betta0 / (alfa0 * alfa0)) <= -h / 2.0)
- {
- IH = I5(-h / 2.0, h / 2.0, delta0 + alfa0, epsilon + (betta0 / alfa0));
- }
- }
- epsilon = epsilonJ(-H, Q, betta);
- betta0 = betta0J(delta, -H, P, betta);
- hi0 = hi0J(-H, Q, r, betta);
- hi1 = sqrt(hi0 - (betta0 * betta0 / (alfa0 * alfa0)));
- epsilon1 = epsilon - (delta0 * betta0 / (alfa0 * alfa0));
- delta1 = delta0 * hi1 / alfa0;
- deltaplus = hi1 + delta1;
- deltaminus = hi1 - delta1;
- dzettaplus = ((alfa0 * h / 2.0) + (betta0 / alfa0)) / hi1;
- dzettaminus = ((-alfa0 * h / 2.0) + (betta0 / alfa0)) / hi1;
- zplus = dzettaplus + sqrt(dzettaplus * dzettaplus + 1.0);
- zminus = dzettaminus + sqrt(dzettaminus * dzettaminus + 1.0);
- if (hi1 > pow(10, -10))
- {
- 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;
- }
- else if (fabs(hi1) <= pow(10, -10))
- {
- if (-h / 2.0 < (-betta0 / (alfa0 * alfa0)) && (-betta0 / (alfa0 * alfa0)) < h / 2.0)
- {
- 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));
- }
- else if ((-betta0 / (alfa0 * alfa0)) >= h / 2.0)
- {
- IminusH = I5(-h / 2.0, h / 2.0, delta0 - alfa0, epsilon - (betta0 / alfa0));
- }
- else if ((-betta0 / (alfa0 * alfa0)) <= -h / 2.0)
- {
- IminusH = I5(-h / 2.0, h / 2.0, delta0 + alfa0, epsilon + (betta0 / alfa0));
- }
- }
- tetta = 1.0 * (IH - IminusH) / betta;
- S0 = S0 + ((etta * tetta * exp(z * xminusy)) / (4 * pi));
- }
- }
- if (abs(x) < 1)
- {
- V0 = exp(z) * sin(abs(x)) / abs(x);
- }
- else if (abs(x) > 1)
- {
- V0 = exp(z * abs(x)) * sin(1) / abs(x);
- }
- //cout << S0 << endl;
- //cout << V0 << endl;
- complex <double> maximum = max(abs(p), abs(S0 - V0) * 1.0 / abs(V0));
- S0 = 0;
- p = maximum;
- }
- }
- cout << "Относительная погрешность=" << abs(p) << endl;
- cout << "______________________" << endl;
- }
- return 0;
- }
Add Comment
Please, Sign In to add comment