Метод верхней релаксации для уравнения Пуассона

[Чекист]

Всем привет! Есть краевая задача для уравнения Пуассона в области, составленной из прямоугольников.

1.jpg
2.jpg

Функции psi_1,3,5(y) = sin(Pi*y), psi_2,4,6(x) = sin(Pi*x), f(x,y) = 1-x*x-y*y.
Параметры a = b = 2, c = d = 1.
Я так понимаю, область можно разбить на 2 прямоугольника и считать для каждого отдельно.
Т.е. получается матрица заполняется граничными условиями вот так?:

3.jpg

Нигде не могу найти нормального описания метода верхней релаксации для данной задачи. Кто знает подскажите пожалуйста. Или ткните носом в гугль.

#define _USE_MATH_DEFINES
#include <iostream>
 
 
int const N = 8;
int const M = 8;
 
                                                            //          |
                                                            //          |
//прототипы                                                 //  ^      b|________
double psi_1(double);                                       //  |i      |        |
double psi_2(double);                                       //  |       |   I    |
double psi_3(double);                                       //          |        |
double psi_4(double);                                       //         d|_ _ _ _ |___________
double psi_5(double);                                       //          |        |           |
double psi_6(double);                                       //          |        |      II   |
double f(double);                                           //          |________|___________|________
                                                            //  ----->           c           a
                                                            //    j
double psi_1(double y){
    return sin(M_PI*y);
}
double psi_2(double x){
    return sin(M_PI*x);
}
double psi_3(double y){
    return sin(M_PI*y);
}
double psi_4(double x){
    return sin(M_PI*x);
}
double psi_5(double y){
    return sin(M_PI*y);
}
double psi_6(double x){
    return sin(M_PI*x);
}
double f(double x, double y){
    return (1 - x*x - y*y);
}
double x[N+1], y[M+1], U[N+1][M+1], u[N+1][M+1];
 
void PuassonPr(void)
{
    //Границы
    int n = N/2;
    int m = M/2;
    
    //I
    for(int j = 0; j <= n ; j++)
    {
        U[0][j] = psi_6(x[j]);
        U[M][j] = psi_4(x[j]);
    }
    for (int i = 0; i <= M; i++)
    {
        U[i][0] = psi_5(y[i]);
        U[i][n] = psi_3(y[i]);
    }
 
    //II
    for(int j = n; j <= N ; j++)
    {
        U[0][j] = psi_6(x[j]);
        U[m][j] = psi_2(x[j]);
    }
    for (int i = 0; i <= m; i++)
    {
        U[i][m] = psi_3(y[i]);
        U[i][N] = psi_3(y[i]);
    }
 
    for (int i = 0; i <= M; i++)
    {
        for (int j = 0; j <= N; j++)
        {
            printf("%.3f \t", U[i][j]);
        }
        std::cout << std::endl;
    }
}
 
int main(int argc, char* argv[]){
    double a = 2.0, d = 1.0;
    double b = 2.0, c = 1.0;
    double Hx = a / N;
    double Hy = b / M;
 
    std::cout.precision(3);
    
    for (int i = 0; i < N; i++)
    {
        x[i] = i*Hx;
    }
    
    for (int j = 0; j < M; j++)
    {
        y[j] = j*Hy;
    }
 
    PuassonPr();
    
    std::system("pause");
    return 0;
}