Для начала — немного старческого брюзжания (хотя нет, получилось много
).
Чтобы правильно решить эту задачу, нужно сочетание сразу трёх вещей: глубокого понимания метода построения разностной схемы, глубокого понимания метода решения системы линейных уравнений, получившихся в результате применения разностной схемы, и, наконец, навыков программирования на выбранном языке (языках). С третьим, вероятно, у вас всё хорошо, но без первых двух этого недостаточно. Глубокое понимание — это не просто понимание получившихся формул, но и понимание, как они были получены, и как можно их проверить, когда вы где-то ошибётесь (а в такой сложной задаче ошибутся все, обязательно и не по одному разу).
При наличии такого понимания вы сможете в отладке посмотреть промежуточные значения и выяснить, в каком месте они начинают отличаться от тех значений, которые должны получиться.
Полностью с вами согласен, что нужно понимание самой задачи. Я изучил теорию данного метода в нескольких учебниках по численным методам, но к сожалению в них изложена информация для тех, кто знаком с теорией и понятием разностных схем (я таковым не являюсь, хоть и был в универе этот предмет). Информация, изложенная в учебнике В. П. Ильина "Методы конечных разностей и конечных объемов для эллиптических уравнений", оказалась для меня немного понятней, хотя метод, к сожалению, я так до конца и не понял.
Вот фрагмент из учебника (стр. 295), использовал формулы (6. 133):
Вот вы пишете, что удалили прогонку по строкам, и результаты стали гораздо лучше. Для меня это удивительно. Вы выводите значения из массива C2, а без выкинутого блока в этом массиве должны остаться только нули начального условия. Возможно, в этом случае вы выводите значения промежуточного шага C1. Но тогда всё равно непонятно, как ненулевые значения попали в столбцы, отличные от столбца 25 (в котором источник). Ведь в методе переменных направлений при прогонке по столбцам все ненулевые значения распространяются только по столбцам, а в 24-м и 26-м столбцах сплошные нули, нули же и должны остаться.
Когда я убрал прогонку по строкам, код выглядит так:
#include <iostream>
#include <fstream>
#include <math.h>
#include <cstdlib>
using namespace std;
const int n = 100;
int i00 = n / 4;
int j00 = n / 4;
//Функция Дирака
double diraca(int x)
{
if (x == 0)
return 1.0;
else
return 0.0;
}
double max(double a, double b)
{
if (a > b)
return a;
else
return b;
}
int main()
{
double V, U;
U = 1.0; V = 1.0;
double D = 0.01;
const double h = 100.0 / double(n);
double a1 = D / (h*h) + U / h;
double b1 = D / (h*h) + V / h;
double c1 = D / (h*h);
double d1 = D / (h*h);
double Q1 = 10.0;
int i, j;
double e1 = 2.0*D / (h*h) + U / h + 2.0*D / (h*h) + V / h;
double e2 = a1 + c1;
double e3 = b1 + d1;
double tau = 1.0;
double f1[n + 1][n + 1];
//решение на n-шаге
double C0[n + 1][n + 1];
//решение на (n+1/2)-шаге
double C1[n + 1][n + 1];
//решение на (n+1)-шаге
double C2[n + 1][n + 1];
//прогоночные коэфф.-ты
double P[n + 1], Q[n + 1];
//начальные условия распределения примеси
for (i = 0; i < n + 1; i++)
for (j = 0; j < n + 1; j++) {
f1[i][j] = Q1 * diraca(i - i00) * diraca(j - j00);
C2[i][j] = 0.0;
C1[i][j] = 0.0;
C0[i][j] = 0.0;
}
double error = 1.0;
while (error > 1e-7)
{
for (i = 0; i < n; i++)
for (j = 0; j < n; j++)
//C0[i][j] = C2[i][j];
C0[i][j] = C1[i][j];
double a, b, c, d[n + 1];
//прогонка по столбцам
for (j = 1; j < n; j++)
{
a = -a1; b = -(1 / tau + e2); c = -c1;
for (i = 1; i < n; i++)
d[i] = b1 * C0[i][j - 1] + d1 * C0[i][j + 1] + (1 / tau - e3)*C0[i][j] + f1[i][j];
P[1] = c / b;
Q[1] = -d[1] / b;
for (i = 2; i < n; i++)
{
P[i] = c / (-a * P[i - 1] + b);
Q[i] = (-d[i] + a * Q[i - 1]) / (-a * P[i - 1] + b);
}
C1[n - 1][j] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
for (i = n - 2; i >= 1; i--)
{
C1[i][j] = P[i] * C1[i + 1][j] + Q[i];
}
}
//граничные условия шаг (n+1/2)
for (i = 0; i < n + 1; i++)
C1[n][i] = C0[n - 1][i];
for (i = 0; i < n + 1; i++)
C1[i][n] = C0[i][n - 1];
for (i = 0; i < n + 1; i++)
C1[0][i] = 0.0;
for (i = 0; i < n + 1; i++)
C1[i][0] = 0.0;
/*
//прогонка по строкам
for (j = 1; j < n; j++)
{
a = -b1; b = -(1 / tau + e3); c = -d1;
for (i = 1; i < n; i++)
d[i] = a1 * C1[j][i - 1] + c1 * C1[j][i + 1] - (1 / tau - e2)*C1[j][i] + f1[j][i];
P[1] = c / b;
Q[1] = -d[1] / b;
for (i = 2; i < n; i++)
{
P[i] = c / (-a * P[i - 1] + b);
Q[i] = (-d[i] + a * Q[i - 1]) / (-a * P[i - 1] + b);
}
C2[j][n - 1] = (d[n - 1] - Q[n - 1] * a) / (-b + P[n - 1] * a);
for (i = n - 2; i >= 1; i--)
{
C2[j][i] = P[i] * C2[j][i + 1] + Q[i];
}
}
//граничные условия шаг (n+1)
for (i = 0; i < n + 1; i++)
C2[n][i] = C1[n - 1][i];
for (i = 0; i < n + 1; i++)
C2[i][n] = C1[i][n - 1];
for (i = 0; i < n + 1; i++)
C2[0][i] = 0.0;
for (i = 0; i < n + 1; i++)
C2[i][0] = 0.0;
*/
error = 0.0;
for (i = 1; i < n; i++)
for (j = 1; j < n; j++)
//error = max(error, fabs(C2[i][j] - C0[i][j]));
error = max(error, fabs(C1[i][j] - C0[i][j]));
cout << error << endl;
}
ofstream out;
out.open("result.txt");
for (i = 1; i < n; i++)
for (j = 1; j < n; j++)
{
out << i << " ";
out << j << " ";
//out << C2[i][j] << "\n";
out << C1[i][j] << "\n";
}
system("pause");
return 0;
}[/code]
Могли бы вы подсказать, если знакомы с данным методом, где может быть ошибка. Я проверял - в формулы данные из условия задачи правильно подставлены. Может я не так, что то понял в методе и есть ошибка где-то в реализации?
02.01.2020, 17:25 
![\[x = 0;C = 0,y = 0;C = 0;x = Lx; \frac{{\partial C}}{{\partial x}} = 0;y = Ly;\frac{{\partial C}}{{\partial y}} = 0;\]](https://dxdy-01.korotkov.co.uk/f/c/0/d/c0d0dbf7e3822859103ee61982cfb27182.png "\[x = 0;C = 0,y = 0;C = 0;x = Lx; \frac{{\partial C}}{{\partial x}} = 0;y = Ly;\frac{{\partial C}}{{\partial y}} = 0;\]")
![\[U = V = 1,D = 0.01,Q = 10,Lx = 100,Ly = 100\]](https://dxdy-01.korotkov.co.uk/f/4/9/9/499d76177f3504d087dd0ac9c5300b5782.png "\[U = V = 1,D = 0.01,Q = 10,Lx = 100,Ly = 100\]")
Разностная схема имеет вид: ![\[U\frac{{{C_{i,j}} - {C_{i - 1,j}}}}{h} + V\frac{{{C_{ij}} - {C_{i,j - 1}}}}{h} = D\left( {\frac{{{C_{i + 1,j}} - 2{C_{i,j}} + {C_{i - 1,j}}}}{{{h^2}}} + \frac{{{C_{i,j + 1}} - 2{C_{i,j}} + {C_{i,j - 1}}}}{{{h^2}}}} \right) + {Q_{i0,j0}}\]](https://dxdy-01.korotkov.co.uk/f/4/8/e/48eda951e406bca58b5a8bb735b9e84782.png "\[U\frac{{{C_{i,j}} - {C_{i - 1,j}}}}{h} + V\frac{{{C_{ij}} - {C_{i,j - 1}}}}{h} = D\left( {\frac{{{C_{i + 1,j}} - 2{C_{i,j}} + {C_{i - 1,j}}}}{{{h^2}}} + \frac{{{C_{i,j + 1}} - 2{C_{i,j}} + {C_{i,j - 1}}}}{{{h^2}}}} \right) + {Q_{i0,j0}}\]")
Нужно решить систему разностных уравнений вида: ![\[- a1{C_{i - 1,j}} - c1{C_{i + 1,j}} - b1{C_{i,j - 1}} - d1{C_{i,j + 1}} + e1{C_{i,j}} = f1,\]](https://dxdy-01.korotkov.co.uk/f/8/4/7/8470d1279658ab166526a5a6de08ac2782.png "\[- a1{C_{i - 1,j}} - c1{C_{i + 1,j}} - b1{C_{i,j - 1}} - d1{C_{i,j + 1}} + e1{C_{i,j}} = f1,\]")
где ![\[a1 = \frac{U}{h} + \frac{D}{{{h^2}}},b1 = \frac{V}{h} + \frac{D}{{{h^2}}},c1 = d1 = \frac{D}{{{h^2}}},\[e1 = \frac{{4D}}{{{h^2}}} + \frac{{U + V}}{h},f1 = {Q_{i0,j0}}\]](https://dxdy-04.korotkov.co.uk/f/b/5/3/b53989f4bcd38bbce6b7b45cfa60ce9382.png "\[a1 = \frac{U}{h} + \frac{D}{{{h^2}}},b1 = \frac{V}{h} + \frac{D}{{{h^2}}},c1 = d1 = \frac{D}{{{h^2}}},\[e1 = \frac{{4D}}{{{h^2}}} + \frac{{U + V}}{h},f1 = {Q_{i0,j0}}\]")
методом переменных направлений Писмана - Речфорда, который описан в учебнике В. П. Ильина "Методы конечных разностей и конечных объемов для эллиптических уравнений" стр. 295. ![\[{C_{0,j}} = 0,j = \overline {0,N} ;{C_{i,0}} = 0,i = \overline {0,N} ;{C_{N,j}} = {C_{N - 1,j}},j = \overline {0,N} ;{C_{i,N}} = {C_{i,N - 1}},i = \overline {0,N} ;\]](https://dxdy-04.korotkov.co.uk/f/3/c/f/3cff23e00488f7d3d2efc7e7d93cf22f82.png "\[{C_{0,j}} = 0,j = \overline {0,N} ;{C_{i,0}} = 0,i = \overline {0,N} ;{C_{N,j}} = {C_{N - 1,j}},j = \overline {0,N} ;{C_{i,N}} = {C_{i,N - 1}},i = \overline {0,N} ;\]")
![\[- a1C_{i - 1,j}^{n - \frac{1}{2}} + \left( {\frac{1}{\tau } + a1 + c1} \right)C_{i,j}^{n - \frac{1}{2}} - c1C_{i + 1,j}^{n - \frac{1}{2}} = b1C_{i,j - 1}^{n - 1} + d1C_{i,j + 1}^{n - 1} + \left( {\frac{1}{\tau } + b1 + d1} \right)u_{i,j}^{n - 1} + {Q_{i0,j0}},\]](https://dxdy-01.korotkov.co.uk/f/0/8/1/08196f3fc079937c3f2ab180114e52db82.png "\[- a1C_{i - 1,j}^{n - \frac{1}{2}} + \left( {\frac{1}{\tau } + a1 + c1} \right)C_{i,j}^{n - \frac{1}{2}} - c1C_{i + 1,j}^{n - \frac{1}{2}} = b1C_{i,j - 1}^{n - 1} + d1C_{i,j + 1}^{n - 1} + \left( {\frac{1}{\tau } + b1 + d1} \right)u_{i,j}^{n - 1} + {Q_{i0,j0}},\]")
а по строкам : ![\[- b1C_{i,j - 1}^n + \left( {\frac{1}{\tau } + b1 + d1} \right)C_{i,j}^n - d1C_{i,j + 1}^n = a1C_{i - 1,j}^{n - \frac{1}{2}} + c1C_{i + 1,j}^{n - \frac{1}{2}} - \left( {\frac{1}{\tau } - a1 - c1} \right)C_{i,j}^{n - \frac{1}{2}} + {Q_{i0,j0}}\]](https://dxdy-03.korotkov.co.uk/f/e/a/5/ea528a04f91b32c37a59398f4479401382.png "\[- b1C_{i,j - 1}^n + \left( {\frac{1}{\tau } + b1 + d1} \right)C_{i,j}^n - d1C_{i,j + 1}^n = a1C_{i - 1,j}^{n - \frac{1}{2}} + c1C_{i + 1,j}^{n - \frac{1}{2}} - \left( {\frac{1}{\tau } - a1 - c1} \right)C_{i,j}^{n - \frac{1}{2}} + {Q_{i0,j0}}\]")
Здесь 
- известное решение, 
, 
- решение на первом полушаге, 
- решение на втором полушаге.
