Дискуссионный математический форумМатематический форум
Математический форум Math Help Planet

Обсуждение и решение задач по математике, физике, химии, экономике

Теоретический раздел
Часовой пояс: UTC + 3 часа [ Летнее время ]
новый онлайн-сервис
число, сумма и дата прописью

Часовой пояс: UTC + 3 часа [ Летнее время ]




Начать новую тему Ответить на тему  [ Сообщений: 4 ] 
Автор Сообщение
 Заголовок сообщения: Решить систему методом переменных направлений
СообщениеДобавлено: 30 дек 2019, 17:26 
Не в сети
Начинающий
Зарегистрирован:
30 дек 2019, 17:06
Сообщений: 3
Cпасибо сказано: 0
Спасибо получено:
0 раз в 0 сообщении
Очков репутации: 1

Добавить очки репутацииУменьшить очки репутации
Здравствуйте! Дано уравнение переноса примеси [math]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;[/math] Граничные условия: [math]x = 0^C = 0,y = 0^C = 0;x = Lx^\frac{{\partial C}}{{\partial x}} = 0;y = Ly^\frac{{\partial C}}{{\partial y}} = 0;[/math] Константы: [math]U = V = 1,D = 0.01,Q = 10,Lx = 100,Ly = 100[/math]. Разностная схема имеет вид: [math]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}}[/math]. Нужно решить систему разностных уравнений вида: [math]- 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[/math], где [math]a1 = \frac{U}{h} + \frac{D}{{{h^2}}},b1 = \frac{V}{h} + \frac{D}{{{h^2}}},c1 = d1 = \frac{D}{{{h^2}}};[/math] методом переменных направлений Писмана - Речфорда, который описан в учебнике В. П. Ильина "Методы конечных разностей и конечных объемов для эллиптических уравнений" стр. 295. Начальные условия: [math]{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} ;[/math]
Расчетные формулы для прогонки по столбцам имеют вид: [math]- 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}},[/math], а по строкам : [math]- 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}}[/math]. Здесь [math]C_{i,j}^{n - 1}[/math] - известное решение, [math]\tau = 1[/math], [math]C_{i,j}^{n - \frac{1}{2}}[/math] - решение на первом полушаге, [math]C_{i,j}^n[/math] - решение на втором полушаге.
Я написал программу для этого метода на 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:
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 градусов):
Изображение

Что не так в программе понять не могу. Прошу вашей помощи.
Обнаружил, что если убрать блок кода с прогонкой по строкам:
#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;
}

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

Почему так происходит не понимаю.

Вернуться к началу
 Профиль  
Cпасибо сказано 
 Заголовок сообщения: Re: Решить систему методом переменных направлений
СообщениеДобавлено: 30 дек 2019, 18:11 
Не в сети
Light & Truth
Аватара пользователя
Зарегистрирован:
15 мар 2016, 15:08
Сообщений: 6677
Cпасибо сказано: 82
Спасибо получено:
1115 раз в 1055 сообщениях
Очков репутации: 187

Добавить очки репутацииУменьшить очки репутации
drovosek писал(а):
Что не так в программе понять не могу.

Кусок вашей программы решает систему линейных уравнений. Он работает? Как проверить? Просто подставьте полученное решение в систему. Должна получиться правая часть.

Вернуться к началу
 Профиль  
Cпасибо сказано 
 Заголовок сообщения: Re: Решить систему методом переменных направлений
СообщениеДобавлено: 04 янв 2020, 14:30 
Не в сети
Начинающий
Зарегистрирован:
30 дек 2019, 17:06
Сообщений: 3
Cпасибо сказано: 0
Спасибо получено:
0 раз в 0 сообщении
Очков репутации: 1

Добавить очки репутацииУменьшить очки репутации
Я подставил данные в вычислительные формулы, так как они записаны в учебнике. Может я неправильно реализовал программу для этого метода? Смотрел теорию по этому методу в разных учебниках, но в них информация изложена сложновато для человека, незнакого с разностными схемами, коим я являюсь, а времени на изучение данного предмета, к сожалению, нет. Программу делал по формулам (6.133) из книги Ильина:
Изображение
Поэтому прошу указать в чем ошибка реализации.

Вернуться к началу
 Профиль  
Cпасибо сказано 
 Заголовок сообщения: Re: Решить систему методом переменных направлений
СообщениеДобавлено: 06 янв 2020, 11:10 
Не в сети
Начинающий
Зарегистрирован:
30 дек 2019, 17:06
Сообщений: 3
Cпасибо сказано: 0
Спасибо получено:
0 раз в 0 сообщении
Очков репутации: 1

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

Вернуться к началу
 Профиль  
Cпасибо сказано 
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему      Страница 1 из 1 [ Сообщений: 4 ]

 Похожие темы   Автор   Ответы   Просмотры   Последнее сообщение 
Не сходится метод переменных направлений (ADI)

в форуме Численные методы

proggamer

2

288

08 июл 2016, 14:04

Решить Систему ДУ методом Адамса

в форуме Численные методы

_Vesna_

4

773

27 апр 2011, 19:14

Решить систему методом Гаусса

в форуме Линейная и Абстрактная алгебра

Silas

1

278

23 дек 2011, 22:34

Решить систему ДУ методом исключения

в форуме Дифференциальные и Интегральные уравнения

VinLo

1

327

10 май 2011, 14:47

Решить систему операционным методом

в форуме Комплексный анализ и Операционное исчисление

ful317

0

242

26 май 2013, 16:02

Решить систему д.у. методом исключения

в форуме Дифференциальное исчисление

neverlucky

2

80

13 апр 2020, 22:08

Методом Гаусса решить систему

в форуме Линейная и Абстрактная алгебра

Neelch

2

372

28 ноя 2013, 22:57

Решить систему методом Гаусса

в форуме Линейная и Абстрактная алгебра

stepanka

5

519

17 ноя 2012, 13:49

Решить систему методом Крамера

в форуме Линейная и Абстрактная алгебра

olga_helga

14

174

20 мар 2020, 16:34

Решить систему методом Крамера

в форуме Линейная и Абстрактная алгебра

SonicTheHedgenog

29

1248

22 янв 2015, 06:57


Часовой пояс: UTC + 3 часа [ Летнее время ]



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и гости: 2


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Перейти:  

Яндекс.Метрика

Copyright © 2010-2020 MathHelpPlanet.com. All rights reserved