trud

#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
int nc = 25;
int width = 8, height = 6, r = 3; // Параметры сетки
double h = 0.25; // Шаг
double *gauss(double **matr, double *y, int n)
{
 double d, s;
 double *x = new double [n], *b = new double [n];
 double **a = new double *[n];
 for (int i = 0; i < n; i++)
 a[i] = new double [n];
 for (int i = 0; i < n; i++)
 for (int j = 0; j < n; j++)
 a[i][j] = matr[i][j];
 for (int i = 0; i < n; i++)
 b[i] = y[i];

 for (int k = 1; k < n; k++) //прямой ход
 for (int j = k + 1; j < n; j++)
 {
 d = a[j][k] / a[k][k];
 for (int i = k; i < n; i++)
 a[j][i] = a[j][i] - d * a[k][i];
 b[j] = b[j] - d * b[k];
 }

 for (int k = n - 1; k >= 1; k--) // обратный ход
 {
 d = 0;
 for (int j = k + 1; j < n; j++)
 {
 s = a[k][j] * x[j];
 d = d + s;
 }
 x[k] = (b[k] - d) / a[k][k];
 }
 delete b;
 for (int i = 0; i < n; i++)
 delete [] a[i];
 delete [] a;
 return x;
}
6int main() {
 int ht_h = height / h + 1, wh_h = width / h + 1, r_h = r / h +
1;
 int n = wh_h * ht_h; // количество уравнений
 int i_0 = r_h, j_0 = wh_h - 1 - r_h;
 double *y = new double [n], *x = new double [n];
 double **a = new double *[n];
 for (int i = 0; i < n; i++) a[i] = new double [n];
 FILE *gp = popen("gnuplot -persist", "w");
 if (gp == NULL)
 {
 printf("Ошибка открытия трубы gnuplot\n");
 return 3;
 }

 // ГУ
 for (int i = 0; i < ht_h; i++)
 for (int j = 0; j < wh_h; j++)
 {
 int y_coor = i_0 - i;
 int x_coor = j - j_0;

 if (j == 0 && i == ht_h/2) { //гу 2 рода
 a[i * wh_h + j][i * wh_h + j] = 1 / h;
 a[i * wh_h + j][i * wh_h + j + 1] = -1 /
h;
 y[i * wh_h + j] = 100;
 }

 else if (j == 0) // левая
 {
 a[i * wh_h + j][i * wh_h + j] = 1;
 y[i * wh_h + j] = 0;
 } else if (i == ht_h - 1) //нижняя
 {
 a[i * wh_h + j][i * wh_h + j] = -1 / h - 1;
 a[i * wh_h + j][(i - 1) * wh_h + j] = 1 / h;
 y[i * wh_h + j] = 0;

 }
 else if ( i < i_0 && j > j_0 && //обработка границ
неправильной формы
 ((pow(y_coor + 1, 2) + pow(x_coor, 2) >= pow(r_h, 2)
||
 pow(y_coor, 2) + pow(x_coor + 1, 2) >= pow(r_h, 2))
&&
 pow(y_coor, 2) + pow(x_coor, 2) < pow(r_h, 2)))
 {
 double l, m;
 if ((l = (sqrt(r_h * r_h - pow(x_coor, 2)) -
y_coor)) > 1) l = 1;
7 if ((m = (sqrt(r_h * r_h - pow(y_coor, 2)) -
x_coor)) > 1) m = 1;

 a[i * wh_h + j][i * wh_h + j] = -1 /(m * h * h) -
1 / (1 * h * h) - 1;
 a[i * wh_h + j][i * wh_h + j + 1] = 1 / (m * (m +
1) * h * h);
 a[i * wh_h + j][i * wh_h + j - 1] = 1 / ((m + 1)
* h * h);
 a[i * wh_h + j][(i + 1) * wh_h + j] = 1 / ((l +
1) * h * h);
 a[i * wh_h + j][(i - 1) * wh_h + j] = 1 / (l *(l
+ 1) * h * h);
 }
 else if (i < i_0 && j > j_0 && //округлая
 ((pow(y_coor - 1, 2) + pow(x_coor, 2) < pow(r_h,
2) ||
 pow(y_coor, 2) + pow(x_coor - 1, 2) < pow(r_h,
2)) &&
 pow(y_coor, 2) + pow(x_coor, 2) >= pow(r_h, 2)))
 {
 a[i * wh_h + j][i * wh_h + j] = 1;
 y[i * wh_h + j] = 200;
 }
 else if (i < i_0 && j > j_0 && //пустая часть
 (pow(y_coor - 1, 2) + pow(x_coor, 2) >=
pow(r_h, 2) &&
 pow(y_coor, 2) + pow(x_coor - 1, 2) >=
pow(r_h, 2) &&
 pow(y_coor, 2) + pow(x_coor, 2) > pow(r,
2)))
 {
 a[i * wh_h + j][i * wh_h + j] = 1;
 y[i * wh_h + j] = 0;
 }
 else if (j == wh_h - 1) //правая
 {
 a[i * wh_h + j][i * wh_h + j] = 1;
 y[i * wh_h + j] = 200;
 }

 else if (i == 0) // верхняя
 {
 a[i * wh_h + j][i * wh_h + j] = 1;
 y[i * wh_h + j] = 200;
 }
 else //остальные узлы
 {
 a[i * wh_h + j][i * wh_h + j] = -2 / (h * h) - 2 /
(h * h) - 1;
 a[i * wh_h + j][i * wh_h + j + 1] = 1 / (h * h);
 a[i * wh_h + j][i * wh_h + j - 1] = 1 / (h * h);
8 a[i * wh_h + j][(i + 1) * wh_h + j] = 1 / (h * h);
 a[i * wh_h + j][(i - 1) * wh_h + j] = 1 / (h * h);
 }
 }

 for (int k = 1; k <= nc; k++)
 {
 for (int i = 0; i < ht_h; i++)
 for (int j = 0; j < wh_h; j++)
 {
 int y_coor = i_0 - i;
 int x_coor = j - j_0;

 if (j > 0 && i > 0 && j < wh_h - 1 && i < ht_h - 1
&&
 (i >= i_0 || j <= j_0))
 y[i * wh_h + j] = -x[i * wh_h + j];
 else if (i < i_0 && j > j_0 &&
 pow(y_coor, 2) + pow(x_coor, 2) <
pow(r_h, 2))
 y[i * wh_h + j] = -x[i * wh_h + j];
 }
 x = gauss(a, y, n);
 }

 // Вывод
 for (int i = 0; i < ht_h; i++)
 {
 for (int j = 0; j < wh_h; j++)
 cout << x[i * wh_h + j] << " ";
 cout << endl;
 }
 fprintf(gp, "set pm3d map\n»);
 fprintf(gp, "set autoscale fix\n");
 fprintf(gp, "set cbtics 25\n");
 fprintf(gp, "set xlabel \"%d\" textcolor lt 6\n", width);
 fprintf(gp, "set label \"%d\" at graph -0.11, graph 0.47
textcolor lt 6\n", height);
 fprintf(gp, "set label \"%d\" at graph 0.9, graph 1.07
textcolor lt 6\n", r);
 fprintf(gp, "splot '-' matrix\n");
 for (int i = ht_h - 1; i >= 0; i--)
 {
 for (int j = 0; j < wh_h; j++)
 {
 fprintf(gp, "%lf ", x[i * wh_h + j]);
 }
 fprintf(gp, "\n");
 }
 fprintf(gp, "e\n");
 fclose(gp);
 return 0;
9}