Решить систему методом переменных направлений

[artey]

Здравствуйте, всех С Новым Годом! Задача такова: дано уравнение переноса примеси $U\frac{{\partial C}}{{\partial x}} + V\frac{{\partial C}}{{\partial y}} = D\left( {\frac{{{\partial ^2}C}}{{\partial {x^2}}} + V\frac{{{\partial ^2}C}}{{\partial {y^2}}}} \right) + Q\delta (x - {x_0})\delta (y - {y_0});0 < x < Lx,$ $0 < y < Ly$
Граничные условия:

\[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\]

Разностная схема имеет вид:

\[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,\]

где

\[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} ;\]

Расчетные формулы для прогонки по столбцам имеют вид:

\[- 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}}\]

Здесь $C_{i,j}^{n - 1}$ - известное решение, $\tau = 1$, $C_{i,j}^{n - \frac{1}{2}}$ - решение на первом полушаге, $C_{i,j}^n$ - решение на втором полушаге.
Я написал программу для этого метода на C++:

Используется синтаксис C++

#include <iostream>
#include <fstream>
#include <math.h>
#include <cstdlib>
using namespace std;
const int n = 100;
//координаты источника
int i00 = n / 4;
int j00 = n / 4;
//Функция Дирака; 1 - там, где источник
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];
 
        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]));
        
 
        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";
        }
    system("pause");
 
    return 0;
}

Построил график решения (распространения примеси) по данным из файла result.txt с помощью скрипта на Python:

Используется синтаксис Python

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
 
path = ".../result.txt"
cols = [0,1,2]
with open(path) as file:
    coord = np.loadtxt(file,usecols=cols)
    file.close()
x = coord[:, 0]
y = coord[:, 1]
z = coord[:, 2]
minv = min(z)
maxv = max(z)
print('min=',minv)
print('max=',maxv)
xi = yi = coord[0:100, 1]
#каждой точке сетки присвоить соответствующее значение
zi = griddata((x, y), z,(xi[None,:], yi[:,None]))
fig, ax = plt.subplots()
cs = ax.contourf(xi,yi,zi,cmap=plt.cm.jet)
cbar = fig.colorbar(cs,label = 'Values')
plt.xlim(min(xi),max(xi))
plt.ylim(min(xi),max(xi))
ax.set_xlabel('X')
ax.set_ylabel('Y')
plt.show()

Вот такой график получился:

На самом деле должно получится что-то вроде этого (так сказал преподаватель), то есть распространение примеси должно быть под углом 45 градусов:

Где ошибка в программе не понимаю.
Обнаружил только, что если убрать блок кода с прогонкой по строкам:

Используется синтаксис C++

//прогонка по строкам
        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;

то получается вот такой вот результат (что более-менее похоже на тот рисунок, где угол 45):

Но это же неправильно, так как противоречит теории для этого метода.

[worm2]

Для начала — немного старческого брюзжания (хотя нет, получилось много ).
Чтобы правильно решить эту задачу, нужно сочетание сразу трёх вещей: глубокого понимания метода построения разностной схемы, глубокого понимания метода решения системы линейных уравнений, получившихся в результате применения разностной схемы, и, наконец, навыков программирования на выбранном языке (языках). С третьим, вероятно, у вас всё хорошо, но без первых двух этого недостаточно. Глубокое понимание — это не просто понимание получившихся формул, но и понимание, как они были получены, и как можно их проверить, когда вы где-то ошибётесь (а в такой сложной задаче ошибутся все, обязательно и не по одному разу).
При наличии такого понимания вы сможете в отладке посмотреть промежуточные значения и выяснить, в каком месте они начинают отличаться от тех значений, которые должны получиться.
Вот вы пишете, что удалили прогонку по строкам, и результаты стали гораздо лучше. Для меня это удивительно. Вы выводите значения из массива C2, а без выкинутого блока в этом массиве должны остаться только нули начального условия. Возможно, в этом случае вы выводите значения промежуточного шага C1. Но тогда всё равно непонятно, как ненулевые значения попали в столбцы, отличные от столбца 25 (в котором источник). Ведь в методе переменных направлений при прогонке по столбцам все ненулевые значения распространяются только по столбцам, а в 24-м и 26-м столбцах сплошные нули, нули же и должны остаться.

[Pphantom]

А я добавлю модераторское брюзжание: для сколько-нибудь длинного кода лучше использовать тэг подсветки синтаксиса, иначе результат почти нечитаем. Выше поправлено.

[artey]

worm2 в сообщении #1433210 писал(а):

Для начала — немного старческого брюзжания (хотя нет, получилось много ).
Чтобы правильно решить эту задачу, нужно сочетание сразу трёх вещей: глубокого понимания метода построения разностной схемы, глубокого понимания метода решения системы линейных уравнений, получившихся в результате применения разностной схемы, и, наконец, навыков программирования на выбранном языке (языках). С третьим, вероятно, у вас всё хорошо, но без первых двух этого недостаточно. Глубокое понимание — это не просто понимание получившихся формул, но и понимание, как они были получены, и как можно их проверить, когда вы где-то ошибётесь (а в такой сложной задаче ошибутся все, обязательно и не по одному разу).
При наличии такого понимания вы сможете в отладке посмотреть промежуточные значения и выяснить, в каком месте они начинают отличаться от тех значений, которые должны получиться.

Полностью с вами согласен, что нужно понимание самой задачи. Я изучил теорию данного метода в нескольких учебниках по численным методам, но к сожалению в них изложена информация для тех, кто знаком с теорией и понятием разностных схем (я таковым не являюсь, хоть и был в универе этот предмет). Информация, изложенная в учебнике В. П. Ильина "Методы конечных разностей и конечных объемов для эллиптических уравнений", оказалась для меня немного понятней, хотя метод, к сожалению, я так до конца и не понял.
Вот фрагмент из учебника (стр. 295), использовал формулы (6. 133):

worm2 в сообщении #1433210 писал(а):

Вот вы пишете, что удалили прогонку по строкам, и результаты стали гораздо лучше. Для меня это удивительно. Вы выводите значения из массива C2, а без выкинутого блока в этом массиве должны остаться только нули начального условия. Возможно, в этом случае вы выводите значения промежуточного шага C1. Но тогда всё равно непонятно, как ненулевые значения попали в столбцы, отличные от столбца 25 (в котором источник). Ведь в методе переменных направлений при прогонке по столбцам все ненулевые значения распространяются только по столбцам, а в 24-м и 26-м столбцах сплошные нули, нули же и должны остаться.

Когда я убрал прогонку по строкам, код выглядит так:

Используется синтаксис C++

#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]

Могли бы вы подсказать, если знакомы с данным методом, где может быть ошибка. Я проверял - в формулы данные из условия задачи правильно подставлены. Может я не так, что то понял в методе и есть ошибка где-то в реализации?

[artey]

Неужели никому не встречалась подобная задача?