Алгоритм решения уравнения методом прогонки для 3х многого

[IRIKA]

СРОЧНО!!! Помогите пожалуйста кто сможет)))

Дано дифференциальное уравнение, которое описывает распределение радона в пористой среде. Это уравнение расписано по разностной схеме и приведено к трехдиагональному виду.
Помогите пожалуйста составить алгоритм решения этого уравнения методом прогонки для 3х многого случая=)
Вот всё что расписывала, дальше застряла((

Можно решать простым разностным методом или методом расщепления вместо прогонки!

Я пыталась кое-что сделать в программе.
Думаю можно взять вот эти коэффициенты :
$T_{12} = 3.8235\cdot 24 \cdot60\cdot60 $ пеириод полураспада 222Rn, сек.
$\Lambda = \frac{0.693}{T_{12}}$ постоянная рaспада
$De = \frac{\Lambda }{k^2}$ коэф-т диффузии

Для численного исследования распределения радона построена мат. модель, которая включает следующее уравнение:

$${D}_{e}\Delta A(x,y,z)+V\nabla A(x,y,z)=\Lambda A(x,y,z)+\Lambda {A}_{\propto } $$

(формула1)
с граничными условиями

$$A(0)=0$$

$$A(\propto )={A}_{\propto }$$

де $A(x,y,z)$– объемная активность радона, ${D}_{e}$– эффективный коэффициент диффузии радона, $V$ – скорость переноса радона, $\Lambda$ – постоянная распада радона.

Данное уравнение решалось в декарторой системе координат методом прогонки. Рассмотриваем уравнение в виде:

$${D}_{e}(\frac{{d}^{2}A}{d{x}^{2}}+\frac{{d}^{2}A}{d{y}^{2}}+\frac{{d}^{2}A}{d{z}^{2}})+V(\frac{dA}{dx}+\frac{dA}{dy}+\frac{dA}{dz})+\Lambda {A}_{\propto }- \Lambda (A(x)+A(y)+A(z))=0$$

фиксируем координаты $Y$, $Z$:

$${D}_{e}\frac{{d}^{2}A}{d{x}^{2}}+V\frac{dA}{dx}++\Lambda {A}_{\propto }- \Lambda A(x)=0$$

фиксируем координаты $x$,$ z$:

$${D}_{e}\frac{{d}^{2}A}{d{y}^{2}}+V\frac{dA}{dy}+\Lambda {A}_{\propto }- \Lambda A(y)=0$$

фиксируем координаты $x$, $y$:

$${D}_{e}\frac{{d}^{2}A}{d{z}^{2}}+V\frac{dA}{dz}++\Lambda {A}_{\propto }- \Lambda A(z)=0$$

таким образом, получим уравнение вдоль одной оси

$${D}_{e}\frac{{A}_{i+1,j,k}-2{A}_{i,j,k}+{A}_{i-1,j,k}}{\Delta {x}^{2}}+V\frac{{A}_{i+1,j,k}+{A}_{i-1,j,k}}{2\Delta x}-\Lambda {A}_{i,j,k}=-\Lambda {A}_{\propto }$$

$${D}_{e}\frac{{A}_{i,j+1,k}-2{A}_{i,j,k}+{A}_{i,j-1,k}}{\Delta {y}^{2}}+V\frac{{A}_{i,j+1,k}+{A}_{i,j-1,k}}{2\Delta y}-\Lambda {A}_{i,j,k}=-\Lambda {A}_{\propto }$$

$${D}_{e}\frac{{A}_{i,j,k+1}-2{A}_{i,j,k}+{A}_{i,j,k-1}}{\Delta {z}^{2}}+V\frac{{A}_{i,j,k+1}+{A}_{i,j,k-1}}{2\Delta z}-\Lambda {A}_{i,j,k}=-\Lambda {A}_{\propto }$$

Обозначим через:

$${A}_{i+1,j,k}={y}_{i+1}, {A}_{i,j,k}={y}_{i}, {A}_{i-1,j,k}={y}_{i-1}$$

$${A}_{i}{y}_{i-1}-{C}_{i}{y}_{i}+{B}_{i}{y}_{i+1}=-{F}_{i}$$

$${A}_{i}=\frac{{D}_{e}}{\Delta {x}^{2}}-\frac{V}{2\Delta x}, {B}_{i}=\frac{{D}_{e}}{\Delta {x}^{2}}+\frac{V}{2\Delta x}, {C}_{i}=\frac{{2D}_{e}}{\Delta {x}^{2}}+\Lambda, {F}_{i}=\Lambda {A}_{\propto }$$

$${A}_{j}=\frac{{D}_{e}}{\Delta {y}^{2}}-\frac{V}{2\Delta y}, {B}_{j}=\frac{{D}_{e}}{\Delta {y}^{2}}+\frac{V}{2\Delta y}, {C}_{j}=\frac{{2D}_{e}}{\Delta {y}^{2}}+\Lambda, {F}_{j}=\Lambda {A}_{\propto }$$

аналогично для $k$.

СИСТЕМА:

$${\alpha }_{i+1}=\frac{{B}_{i}}{{C}_{i}-{\alpha }_{i}{A}_i}$$

$${\beta }_{i+1}=\frac{{A}_{i}{\beta }_{i}+{F}_{i}}{{C}_{i}-{\alpha}_{i}{A}_i}$$

$${\alpha }_{1}=0, {\beta }_{1}=0$$

$${y}_{i}={\alpha }_{i+1}{y}_{i+1}+{\beta }_{i+1}, i=N-1,N-2,...0$$

$${y}_{N}={A}_{\propto }$$

[Pphantom]

Тема перемещена из форума «Физика» в форум «Карантин»
по следующим причинам:

- неправильно набраны формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы);
- отсутствуют собственные содержательные попытки решения задачи.

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.

[Lia]

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»

[svv]

Претензии к уравнению.

0. Каждая производная $A$ по координате должна быть частной: $\frac{\partial A}{\partial x}$, в коде пишется \frac{\partial A}{\partial x}.

1. $V\nabla A$ — такого слагаемого быть не может.
У Вас $A$ скалярная величина, концентрация радона (по-моему, более уместный термин, чем «активность»). Тогда $\nabla A$ — это будет градиент концентрации, вектор с компонентами $(\frac{\partial A}{\partial x}, \frac{\partial A}{\partial y}, \frac{\partial A}{\partial z})$. Это имеет смысл. Теперь,
$\bullet$ Если скорость — скаляр $V$, то $V\nabla A$ даёт вектор, что противоречит скалярности остальных слагаемых. Складывать вектор со скалярами бессмысленно.
$\bullet$ Если скорость — вектор $\vec V$ с компонентами $(V_x, V_y, V_z)$, тогда надо писать скалярное произведение: $\vec V\cdot\nabla A$ или $(\vec V, \nabla A)$. Но в таком случае Вы неправильно расписали это ниже, должно быть $V_x\frac{\partial A}{\partial x}+V_y\frac{\partial A}{\partial y}+V_z\frac{\partial A}{\partial z}$.

Выражение $\frac{\partial A}{\partial x}+\frac{\partial A}{\partial y}+\frac{\partial A}{\partial z}$ для скалярной концентрации $A$ никакого физического смысла не имеет.

2. Вы решаете стационарное уравнение, в котором $A$, несмотря на процессы переноса и распада, не зависит от времени в любой точке. Это физически правдоподобно? (Здесь Вы можете просто ответить: «Да, я рассматриваю именно стационарный случай установившейся концентрации»)

3. Хотя бы мысленно, раз нельзя уже подправить, замените «3х многого» на «трёхмерного».

IRIKA в сообщении #1023841 писал(а):

фиксируем координаты $Y$, $Z$:

4. Хорошо, фиксируем. Но от того, что Вы рассматриваете какие-то зависимости от $x$ при фиксированных $y$ и $z$, само уравнение не может измениться. Иначе было бы очень просто решать уравнения с частными производными в многомерных областях. А у Вас с помощью такой фиксации производные по $y$ и $z$ исчезли, и задача стала одномерной. Вы просто выбросили содержательные слагаемые, выражающие зависимость $A$ от $y$ и $z$, потому что они Вам мешали.

Иначе говоря, из полного уравнения никак не следует справедливость трех более простых уравнений с зависимостью $A$ только от $x$, только от $y$ и только от $z$.

IRIKA в сообщении #1023841 писал(а):

Помогите пожалуйста составить алгоритм решения этого уравнения методом прогонки

См. книгу Ильин, Кузнецов. Трехдиагональные матрицы и их приложения.

Итак, Вам нужно выполнить следующие этапы:
$\bullet$ Правильно записать дифференциальное уравнение.
$\bullet$ Правильно построить его разностную аппроксимацию.
$\bullet$ Правильно свести разностное уравнение к трёхдиагональной СЛАУ.
$\bullet$ Составить алгоритм решения трёхдиагональной системы методом прогонки и решить её.

[svv]

Сейчас только заметил.

IRIKA в сообщении #1023841 писал(а):

${D}_{e}(\frac{{d}^{2}A}{d{x}^{2}}+\frac{{d}^{2}A}{d{y}^{2}}+\frac{{d}^{2}A}{d{z}^{2}})+V(\frac{dA}{dx}+\frac{dA}{dy}+\frac{dA}{dz})+\Lambda {A}_{\propto }- \Lambda (A(x)+A(y)+A(z))=0$

Если выражение $\frac{dA}{dx}+\frac{dA}{dy}+\frac{dA}{dz}$ (после исправления полных производных на частные) было только физически бессмысленно, то $A(x)+A(y)+A(z)$ — ещё и математически. Короче говоря, пока не напишете корректное уравнение, о методе прогонки даже не думайте.

[IRIKA]

Всё, что я имела изначально: первое уравнение с граничными условиями.
Нужно вывести решение системы диф. уравнений и сравнить с поверхностями на картинках.








-- 10.06.2015, 13:29 --

Это то, что получилось в прошлом году в МАТЛАБЕ.
На тот момент я использовала точки, в которых измерялся радон и значения, которые мы в этих точках получили.

-- 10.06.2015, 13:56 --

Концентрацию радона меряли мы сами в кабинетах университета. На четвёртой картинке видно 3 этажа и точки в которых измеряли радон. Там где высокая концентрация, кабинет был долго закрыт и радон просто накопился.

[svv]

А откуда Вы брали систему уравнений? (Есть ли ссылочка?)
Кстати, в начале же было одно уравнение, а не система?
Что такое скорость переноса? Я имею в виду, почему недостаточно коэффициента диффузии, есть ещё какой-то механизм переноса? Откуда она известна?

[IRIKA]

http://cyberleninka.ru/article/n/diffuz ... kih-sredah

Вот ссылка статьи, опираясь на те формулы я и создавала эти формулы.

-- 10.06.2015, 15:36 --

мне нужно создать опять три поверхности и сравнить их с теми что на картинках, только тогда всё было проще, я знала и точки и значения в точках, а теперь это всё нужно сделать с использованием уравнения, вот я его и решала. Как я считаю огромная проблема еще в том, что многие коэффициенты мне не известны и я должна сама их подобрать.

[svv]

Понятно, в чём дело — там одномерная задача. И у Вас неверно сделан переход к трехмерному случаю.
В качестве первого приближения (хотя мне ещё надо будет подумать) трёхмерное обобщение должно выглядеть так:
$D\Delta A-\mathbf v\cdot\nabla A-\lambda A+\lambda A_{\infty}=0$
Тут много идеологических вопросов, например, откуда Вам брать поле скоростей $\mathbf v$ — это гораздо сложнее, чем подставить скалярную константу. Скорость конвекции может быть и константой, но — векторной. Я думаю, векторность $\mathbf v$ достаточно очевидна: в принципе, некий существующий в данной точке «сквозняк» может сдувать атомы радона не только с различной скоростью, но и в произвольном направлении.

Судя по графикам, Вам необходима именно трёхмерная модель в задаче?

(исправил знак)

[IRIKA]

Да, именно трёхмерной.

-- 10.06.2015, 15:46 --

Там по началу кажется всё более менее понятно, но когда начинаешь расписывать всё в трёхмерке с каждым шагом всё больше и больше вопросов возникает.

-- 10.06.2015, 15:49 --

Всю информацию про радон я узнавала с методичек и т.д.
Конкретного человека, который бы рассказал про скорости и всё остальное нет.

[svv]

А откуда бы Вы брали скалярное значение скорости переноса?
Для меня это, буквально, средняя скорость сквозняка в помещении.

[IRIKA]

Как и многие коэффициенты, я бы их с потолка.
К сожалению этих значений нет, в этом для меня огромная проблема, большенство значений нужно подобрать самим.

-- 10.06.2015, 16:12 --

брала бы их с потолка

[IRIKA]

Можно решать простым разностным методом или методом расщепления вместо прогонки!

-- 10.06.2015, 22:53 --

Я пыталась кое-что сделать в программе.
Думаю можно взять вот эти коэффициенты :
$T_{12} = 3.8235\cdot 24 \cdot60\cdot60 $ пеириод полураспада 222Rn, сек.
$\Lambda = \frac{0.693}{T_{12}}$ постоянная рaспада
$De = \frac{\Lambda }{k^2}$ коэф-т диффузии

-- 10.06.2015, 23:00 --



Что делать со скоростью $V$ я не знаю((

[vasya321]

Для решения уравнения

$$ \nabla \left( D(z) \nabla A(x,y,z)\right) - \mathbf{V} \, \nabla A(x,y,z) - \lambda ( A(x,y,z) - A_{\infty} ) = 0 $$

советую использовать конечно-объемный метод. Данный метод подробно рассмотрен в книге "Патанкар - Численные методы решения задач теплообмена и динамики жидкости". На стр. 84 приводятся формулы для коэффициент системы лин. алг. уравнений для 3-х мерного случая (лучше, конечно, самому их вывести, чтобы разобраться). Получившуюся СЛАУ можно решить методом Гаусса-Зейделя.

Вообще не ясно:
1) Вам надо рассчитать распределение концентрации радона в доме или в грунте
2) Какие граничные условия на плоскостях $(x,0,z)$, $(x,H_y,z)$, $(0,y,z)$, $(H_x,y,z)$, здесь $H_x$, $H_y$ — координаты границы расчетной области по $x$ и $y$ соответственно

[IRIKA]

Концентрация радона в помещении. В кабинетах университета.

-- 11.06.2015, 01:30 --


Вот точки и значения в них:
1т.: 25, 9, 0.5, 206;
2т.: 8, 3, 0.5, 64;
3т.: 9, 50, 0.5, 62;
4т.: 12, 80, 0.5, 80;
5т.: 8, 8, 3.75, 75;
6т.: 12, 20, 3.75, 209;
7т.: 5, 42, 3.75, 26;
8т.: 3, 65, 3.75, 68
9т.: 14, 14, 7, 106
10т.: 12, 38, 7, 32
11т.: 3, 46, 7, 18
12т.: 12, 86, 7, 37

Вот матрица

$$\begin{bmatrix} 25& 9& 0.5& 206\\ 8& 3& 0.5& 64\\ 9& 50& 0.5& 62\\ 12& 80& 0.5& 80\\ 8& 8& 3.75& 75\\ 12& 20& 3.75& 209\\ 5& 42& 3.75& 26\\ 3& 65& 3.75& 68\\ 14& 14& 7& 106\\ 12& 38& 7& 32\\ 3& 46& 7& 18\\ 12& 86& 7& 37\\ \end{bmatrix}$$

Я брала схему расположения кабинетов, мысленно брала оси $x, y$, $z$ - это высота.
Первое значение это $x$
второе - $y$
третье - $z$
четвёртое - это концентрация радона в кабинете.