Численное решение двумерного нестационарного уравнения тепло

[SteelRend]

Здравствуйте. Необходимо численно решить двумерное нестационарное уравнение теплопроводности с переменными коэффициентами (и правой частью) вида
$\frac{\partial u}{\partial \tau} = \frac{\partial}{\partial x_1} (k_1 \frac{\partial u}{\partial x_1} ) + \frac{\partial}{\partial x_2} (k_2 \frac{\partial u}{\partial x_2} ) + f$

$0 < x_1 < L_1, 0 < x_2 < L_2, 0 < t \leqslant T$

$u=u(x_1,x_2,t) $
$f=f(x_1,x_2,t)$
$k_1=k_1(x_1,x_2)$
$k_2=k_2(x_1,x_2)$

где $u$ -- непрерывная функция, удовлетворяющая условиям

(начальные условия)
$u(x_1,x_2,0)=u_0(x_1,x_2)$

(граничные условия)
$\frac{\partial u (0,x_2,t) }{\partial x_1} = g_1 (0,x_2,t)$
$\frac{\partial u (L_1,x_2,t) }{\partial x_1} = g_2 (L_1,x_2,t)$
$\frac{\partial u (x_1,0,t) }{\partial x_2} = g_3 (x_1,0,t)$
$\frac{\partial u (x_1,L_2,t) }{\partial x_2} = g_4 (x_1,L_2,t)$

Как завещал Самарский(Тихонов А.Н., Самарский А. А. Уравнения математической физики), используем (неявную) разностную схему (сетка равномерная)
$\frac{y^{n+1}_{i,j}-y^{n}_{i,j}}{\tau}=\frac{1}{h_1}({a^{(1)}_{i+1,j} \frac{y^{n+1}_{i+1,j}-y^{n+1}_{i,j}}{h_1} - {a^{(1)}_{i,j} \frac{y^{n+1}_{i,j}-y^{n+1}_{i-1,j}}{h_1})} + \frac{1}{h_2}({a^{(2)}_{i,j+1} \frac{y^{n+1}_{i,j+1}-y^{n+1}_{i,j}}{h_2} - {a^{(2)}_{i,j} \frac{y^{n+1}_{i,j}-y^{n+1}_{i,j-1}}{h_2})} + \phi^{n+1}_{i,j}$
$1 \leqslant i \leqslant N_1-1, 1 \leqslant j \leqslant N_2-1, 1 \leqslant n \leqslant K - 1$

($0, ... , N_1, 0, ... , N_2, 0, ... , K$ -- номера узлов сетки по $x_1,x_2,t$ соотв.)

где
${a^{(1)}_{i,j} = k_1((i+\frac{1}{2}) \cdot h_1,j \cdot h_2)$
${a^{(2)}_{i,j} = k_2(i \cdot h_1,{(j+\frac{1}{2})} \cdot h_2)$

для начальных условий:
$y^0_{i,j} = u_0(i \cdot h_1, j \cdot h_2), n = 0$

Как аппроксимировать граничные?

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

Или же следует использовать другой метод?

[GAA]

1. Для двумерной задачи попробуйте использовать продольно-поперечную схему (по-другому — схему переменных направлений;на рус.: Longitudinal-Transverse scheme, на англ.: Alternating Direction Implicit Scheme).
См., например, Калиткин Н.Н. Численные методы. (Легко нагуглить в Сети и почти во всех библиотеках есть.)
2. Для аппроксимации краевых условий второго рода, видимо, проще всего использовать метод «фиктивных узлов» (pdf).
[В п. 2 («Пример реализации схемы переменных направлений») предлагается для аппроксимации условий второго рода строить сетку так, чтобы находились узлы, симметрично отстоящие от границы на $h/2$. Тогда $( u_1-u_0)/h$ будет аппроксимировать производную со вторым порядком $h$.]

-- Вс 02.04.2017 20:04:19 --

В книге Тихонова и Самарского это §4 дополнения I.

[ewert]

SteelRend в сообщении #1206030 писал(а):

$\frac{y^{n+1}_{i,j}-y^{n}_{i,j}}{\tau}=\frac{1}{h_1}({a^{(1)}_{i+1,j} \frac{y^{n+1}_{i+1,j}-y^{n+1}_{i,j}}{h_1} - {a^{(1)}_{i,j} \frac{y^{n+1}_{i,j}-y^{n+1}_{i-1,j}}{h_1})} + \frac{1}{h_2}({a^{(2)}_{i,j+1} \frac{y^{n+1}_{i,j+1}-y^{n+1}_{i,j}}{h_2} - {a^{(2)}_{i,j} \frac{y^{n+1}_{i,j}-y^{n+1}_{i,j-1}}{h_2})} + \phi^{n+1}_{i,j}$

При каждом фиксированном $j$ добавьте к набору уравнений с разными $i$ по краям по одному уравнению -- на частные производные по иксам. Эти производные должны быть трёхточечными, чтоб не потерять 2-й порядок точности схемы. Затем желательно скомбинировать нулевое уравнение с первым и последнее с предпоследним так, чтобы матрица этой системы (отвечающая фиксированному $j$) оказаласи симметричной (и, следовательно, трёхдиагональной).

Затем аналогичные манипуляции надо проделать по лругой переменной. Получится замкнутая линейная система с симметричной матрицей (саму матрицу, естественно, никто и никогда не выписывает в явном виде).

SteelRend в сообщении #1206030 писал(а):

Учитель говорил, что эту задачу можно решать методом Зейделя

Зайдель -- это, конечно, здорово. Однако подобные задачи куда эффективнее решаютя методом переменных направлений = дробных шагов. Схема усложняется, но совсем ненамного, а вот скорость счёта увеличивается очень сильно.

[ewert]

(Оффтоп)

GAA в сообщении #1206034 писал(а):

схему переменных направлений;на рус.: Longitudinal-Transverse scheme

вряд ли всё же это на рус.; скорее на греч.

[GAA]

(О названии)

«Longitudinal-Transverse scheme» я встречал только в статьях (преимущественно переводах русскоязычных на английский и даже просто в резюме) и переводах на английский книг советских авторов. Например, Alexander A. Samarskii The Theory of Difference Schemes.
Легко нагуглить кучу примеров. В англоязычных я не встречал.
Но особо и не искал. Специально этим вопросом не интересовался. Просто указал ключевые слова для поиска в Сети.

[ewert]

(Оффтоп)

GAA в сообщении #1206062 писал(а):

Просто указал ключевые слова для поиска в Сети.

Но они тупо не найдутся при поиске на рус. Я только лишь об этом, просто позабавило.

[GAA]

(2ewert)

Я как раз таким способом нашёл нужную статью на русском языке. (С аннотацией на типа английском.) Но, думаю, ТС и стандартными учебниками вполне удовлетворится; не та задача. А вот для осмысленного названия процедуры может и будет полезно.

[ewert]

Ну, в добивку.

GAA в сообщении #1206034 писал(а):

Для аппроксимации краевых условий второго рода, видимо, проще всего использовать метод «фиктивных узлов»

Это, видимо, имеется в виду примерно следующее (по ссылке я, ессно, не заходил за леностью).

Зафиксируем разностные уравнения для крайних узлов. Кои некорректны, ибо задействуют минус первые узловые значения слева и соотв. справа.

Но вот теперь есть у нас зато и граничные условия. Которые в форме разностных производных симметричного вида выглядят как $\frac{y_{1}-y_{-1}}{2h}=\text{чего-то там насчёт правых частей гр. усл.}$.
Ну а теперь остаётся вычесть друг из дружки крайние уравнения попарно, чтоб исключить фиктивные узловые значения справа и слева.

[SteelRend]

Спасибо за помощь.

Получается, мне нужно добавить узлы с номерами $i=-1,j=0,...,N_2$ и $i=N_1+1,j=0,...,N_2$ (слева и справа), записать уравнение при $i=0$ как если бы эти узлы были внутренними, выразить из уравнения для $i=0$ (которые теперь стали внутренними) $y^{n+1/2}_{-1,j}$, подставить это в граничные условия (для $i=0$), исключив неизвестные с этим номерами, затем выразить $y^{n+1/2}_{0,1}$ через $y^{n+1/2}_{1,1}$, и посмотреть чему будет равен первый прогоночный коэффициент $\alpha_1$?

Для $i=N_1+1,j=0,...,N_2$ записать уравнение при $i=N_1$ как если бы эти узлы были внутренними, выразить из уравнения для $i=N_1$(которые теперь стали внутренними) $y^{n+1/2}_{N_1+1,j}$, подставить это в граничные условия (для $i=N_1$), исключив неизвестные с этим номерами, и затем найти $y^{n+1/2}_{N_1,1}$?

И аналогично для прогонки по столбцам. Так?

И затем то же самое для слоя $n+1$. Правильно?
Ув. ewert, когда Вы говорили "вычесть друг из дружки крайние уравнения попарно", Вы имели в виду вычесть из уравнения для внутренних узлов при $i=0$ (соответственно при $i=N_1$) уравнение для граничных узлов, содержащее "некорректные" значения при $i=-1$ (соответственно $i=N_1+1$)?

[svv]

SteelRend в сообщении #1206030 писал(а):

Необходимо численно решить двумерное нестационарное уравнение теплопроводности с переменными коэффициентами (и правой частью) вида
$\frac{\partial u}{\partial \tau} = \frac{\partial}{\partial x_1} (k_1 \frac{\partial u}{\partial x_1} ) + \frac{\partial}{\partial x_2} (k_2 \frac{\partial u}{\partial x_2} ) + f$

Уравнение неоднородное, нестационарное, всё переменное... Но даже при этом возьмите $k_1=k_2=k$.

Выбор $k_1\neq k_2$ означает, что коэффициент теплопроводности анизотропен — различен в разных направлениях. Но если так, то, скорее всего, всё ещё хуже:
$\frac{\partial u}{\partial \tau} = \sum\limits_{i,j}\frac{\partial}{\partial x_i} (k_{ij} \frac{\partial u}{\partial x_j} ) + f =$
$=\frac{\partial}{\partial x_1} (k_{11} \frac{\partial u}{\partial x_1} ) +\frac{\partial}{\partial x_1} (k_{12} \frac{\partial u}{\partial x_2} ) +\frac{\partial}{\partial x_2} (k_{21} \frac{\partial u}{\partial x_1} ) +\frac{\partial}{\partial x_2} (k_{22} \frac{\partial u}{\partial x_2} ) +f$
Единственное слабое утешение: $k_{12}=k_{21}$.

[SteelRend]

svv
В $§ 4$ дополнения I из "Уравнений математической физики" Самарского и Тихонова встречается то самое уравнение, которое я выписал в заглавном сообщении. Такое же уравнение я видел в книге Годунова по разностным схемам. На это же уравнение мне указывал учитель. В остальной литературе там один и тот же переменный коэффициент, либо постоянный. И я не знаю, каким "физическим условиям" процесса переноса тепла это уравнение соответствует. Кстати, $k_1, k_2$ -- коэффициенты температуропроводности.

[svv]

OK.

[SteelRend]

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

$y^{n+1/2}_{0,j}=\frac{y^{n}_{0,j}+y^{n+1}_{0,j}}{2} - \frac{\tau}{4}\Lambda_2(y^{n+1}_{0,j}+y^{n}_{0,j})$

однако $y^{n+1}_{0,j}$ в моём случае так легко не определяются, поскольку граничное условие 2-го рода. Что посоветуете? Может, я что-то не правильно понял?

[ewert]

SteelRend в сообщении #1206080 писал(а):

когда Вы говорили "вычесть друг из дружки крайние уравнения попарно", Вы имели в виду вычесть

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

SteelRend в сообщении #1206080 писал(а):

$y^{n+1/2}_{N_1,1}$

Забудьте пока про дробные шаги, разберитесь с вашей неявной схемой.

Вообще же метод дробных шагов -- штука идейно простая. У Вас в правой части стоит сумма двух разностных операторов -- один отвечает второй производной по первой переменной (включая граничные условия), другой -- по второй. И оба применяются к решению на новом временнОм слое. В явной схеме оба применялись бы к решению на слое старом. Ну а в схеме переменных направлений на первом полушаге первая половина оператора применяется к значениям на старом слое, вторая -- к новым. На втором полушаге -- наоборот.

В чём выигрыш? Как и в неявной схеме, на каждом шаге приходится решать системы уравнений. Но, в отличие от неявной схемы, теперь на каждом полушаге получается набор независимых систем по одной переменной при фиксированной другой. Эти системы трёхдиагональны и потому эффективно решаются методом прогонки; общее время счёта $O(n^2)$ на шаг, т.е. минимально возможное (по порядку). При этом устойчивость остаётся абсолютной, как и в неявной схеме. Зато плюс дополнительный бонус: теперь из-за определённой симметричности по времени погрешность оценивается как $O(\tau^2+h^2)$ (в отличие от явной и неявной схем, для которых всего лишь $O(\tau+h^2)$). Поэтому можно брать $\tau\sim h$ вместо $\tau\sim h^2$, что сильно ускоряет работу при той же точности.

svv в сообщении #1206085 писал(а):

Выбор $k_1\neq k_2$ означает, что коэффициент теплопроводности анизотропен — различен в разных направлениях. Но если так, то, скорее всего, всё ещё хуже

А это не имеет значения -- главное, что дифоператор симметричен и знакоопределён (правда, вырожден, поскольку задача Неймана, ну и бог с ней). И это же при грамотной реализации переносится и на оператор разностный.

Правда, если бы анизотропия была бы полной (тензор теплопроводностей не был бы диагонален), то это помешало бы схеме переменных направлений. Но всё остальное сохранилось бы.

-- Пн апр 03, 2017 11:12:45 --

SteelRend в сообщении #1206089 писал(а):

Кстати, $k_1, k_2$ -- коэффициенты температуропроводности.

Кстати, теплопроводности. Температуропроводность получится, если разделить на теплоёмкость (которая в вашей схеме подразумевается единичной). Но это в неоднородной среде несколько бессмысленно.

[SteelRend]

Значит, уравнение для левой границы

$\frac{y^{n+1/2}_{1,j} - y^{n+1/2}_{-1,j} }{2 h_1} = g^{n+1/2}_{0,j} (*)$

Уравнение для $i=0$

$\frac{y^{n+1/2}_{0,j} - y^{n}_{0,j} }{0,5 \tau} = \frac{y^{n+1/2}_{1,j} - 2y^{n+1/2}_{0,j} + y^{n+1/2}_{-1,j} }{h^2_1$} + $\Lambda_2 y^{n}_{0,j}+\phi^{n+1/2}_{0,j} (**)$


Из этих двух уравнений исключить? Выражаю $y^{n+1/2}_{-1,j}$ из $(*)$, подставляю в $(**)$, получаю уравнение для $i=0$, не содержащее значений для фиктивных узлов.

Далее я выражаю $y^{n+1/2}_{0,j}$ для прогонок по $x_1$

$y^{n+1/2}_{0,j}=\frac{2}{h_1} \frac{1}{m} y^{n+1/2}_{1,j}+\frac{1}{m}(-\frac{2}{h_1} g^{n+1/2}_{0,j} + \frac{1}{0,5 \tau} y^{n}_{0,j} + \Lambda_2 y^{n}_{0,j} + \phi^{n+1/2}_{0,j})$
$m=\frac{1}{0,5 \tau}+\frac{2}{h^{2}_1}}$

И тогда

$\alpha_1 = \frac{2}{h_1} \frac{1}{m}$ ,
$\beta_1 = \frac{1}{m}(-\frac{2}{h_1} g^{n+1/2}_{0,j} + \frac{1}{0,5 \tau} y^{n}_{0,j} + \Lambda_2 y^{n}_{0,j} + \phi^{n+1/2}_{0,j})$

Верно?