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

#include "pch.h"
#include <iostream>
#include <conio.h>
#include <string>
#include <cstring>
#include <math.h>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <algorithm>
using namespace std;

int main()
{
    int M = 25, N = 25;
    long double V0, K0 = 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, xminusy, R;
    double A = N * h, B = M * H;
    //double uq = pi / (1.0 * N), vl = pi / (2.0 * M);
    setlocale(LC_ALL, "Rus");
    for (int i = 1; i < 6; i++)
    {
        R = 1 + pow(10, (-1)*i);
        //cout << "............................." << endl;
        cout << "R=" << R << endl;
        for (int q = 0; q <= 2 * N; 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;
                //cout << "x1=" << x1 << endl;
                //cout << "x2=" << x2 << endl;
                //cout << "x3=" << x3 << endl;
                x = pow(pow(x1, 2) + pow(x2, 2) + pow(x3, 2), 0.5);
                //cout << "x=" << x << endl;     //верно

                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);
                        etta = sin(vm);
                        xminusy = pow(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2), 0.5);
                        //cout << "1:" << xminusy << endl;
                        K0 = K0 + ((h * H * etta * 1.0  * cos(un) * sin(vm) / xminusy) / (4 * pi));
                        //xminusy = sqrt(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2));
                        //cout << "2:" << xminusy << endl;
                    }
                }
                //cout << "K0=" << K0 << endl;
                if (abs(x) < 1)
                {
                    V0 = abs(x) * x1 * 1.0 / (3 * R);
                }
                else if (abs(x) > 1)
                {
                    V0 = x1 / (3 * R * 1.0 * pow(abs(x), 2));
                }

                //cout << "V0=" << V0 << endl;
                long double maximum = max(p, abs(K0-V0));
                K0 = 0;
                p = maximum;
            }
        }
        cout << "Абсолютная погрешность=" << p << endl;
        cout << "______________________" << endl;
    }
    return 0;
}





#include "pch.h"
#include <iostream>
#include <conio.h>
#include <string>
#include <cstring>
#include <math.h>
#include <cstdlib>
#include <ctime>
#include <cmath>
using namespace std;

int main()
{
    int M = 100, N = 100;
    long double V0, K0 = 0, S0 = 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, xminusy, R;
    double A = N * h, B = M * H;
    double tetta, betta, y1pr, y2pr, y3pr;
    double uq = pi / (1.0 * N), vl = pi / (2.0 * M);
    for (int i = 1; i < 6; i++)
    {
        R = 1 + pow(10, (-1)*i);
        cout << "R=" << R << endl;
        y11 = sin(vl)*cos(uq);
        y21 = sin(vl)*sin(uq);
        y31 = cos(vl);
        x1 = y11 * R;
        x2 = y21 * R;
        x3 = y31 * R;
        cout << "x1=" << x1 << endl;
        cout << "x2=" << x2 << endl;
        cout << "x3=" << x3 << endl;
        x = pow(pow(x1, 2) + pow(x2, 2) + pow(x3, 2), 0.5);
        cout << "x=" << x << endl;
        // Вычисление K0
        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);
                etta = sin(vm);
                xminusy = pow(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2), 0.5);
                K0 = K0 + ((h * H * etta * 1.0  * cos(un) * sin(vm) / xminusy) / (4 * pi));
            }
        }
        // Вычисление S0
        /* 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;
                y1pr = cos(un) * cos(vm);
                y2pr = sin(un) * cos(vm);
                y3pr = (-1)*sin(vm);
                betta = sqrt(pow(y1pr, 2) + pow(y2pr, 2) + pow(y3pr, 2));

                // Вычисление I(H)






                tetta = 1.0 / betta;
                y1 = sin(vm)*cos(un);
                y2 = sin(vm)*sin(un);
                y3 = cos(vm);
                etta = sin(vm);
                xminusy = pow(pow(x1 - y1, 2) + pow(x2 - y2, 2) + pow(x3 - y3, 2), 0.5);
                S0 = S0 + (etta * tetta);
            }
        } */

        cout << "K0=" << K0 << endl;
        if (abs(x) < 1)
        {
            V0 = abs(x) * x1 * 1.0 / (3 * R);
        }
        else if (abs(x) > 1)
        {
            V0 = x1 / (3 * R * 1.0 * pow(abs(x), 2));
        }

        cout << "V0=" << V0 << endl;
        cout << "______________________" << endl;
    }

    return 0;
}