Численное решение уравнения переноса

[Pete[r]]

Имеется ур-е конвективного переноса $\dfrac{\partial\varphi}{\partial t} + u\dfrac{\partial\varphi}{\partial x} = 0$, где $u$ --- скорость, $\varphi$ --- какая-то скалярная функция.
На сетке $0\leqslant x = ih\leqslant s = Nh$, $0\leqslant t = j\tau\leqslant M\tau$ надо использовать такую четырёхточечную схему: $\dfrac{\varphi_i^{j+1} - \varphi_i^j}{\tau} + u\dfrac{\varphi_{i+1}^{j+1} - \varphi_{i-1}^{j+1}}{2h} = 0$.

$$ \xymatrix{ (i-1, j+1) \ar@{-}[r] & (i, j+1) \ar@{-}[r] \ar@{-}[d] & (i+1, j+1)\\ & (i, j) &\\ } $$

Эта схема сводится к СЛАУ

$$ \begin{cases} A_if_{i-1} - C_if_i + B_if_{i+1} = -F_i,\\ i = 1,\ldots, N, \end{cases} $$

$A_i = -\dfrac{\tau u}{2h}$, $B_i = \dfrac{\tau u}{2h}$, $C_i = -1$, $F_i = \varphi_i^j$, $f_{i-1} = \varphi_{i-1}^{j+1}$, $f_i = \varphi_i^{j+1}$, $f_{i+1} = \varphi_{i+1}^{j+1}$.

Коэффициенты образуют трёхдиагональную матрицу, которая решается методом прогонки следующим образом.
Граничные условия задаём в виде $f_0 = \eta_1f_1 + \mu_1$, $f_N = \eta_2f_{N-1} + \mu_2$.
Получаем алгоритм:

$$ \begin{cases} \alpha_1 = \eta_1, \beta_1 = \mu_1,\\ \alpha_{i+1} = \dfrac{B_i}{C_i - \alpha_iA_i}, \beta_{i+1} = \dfrac{A_i\beta_i + F_i}{C_i - \alpha_iA_i}, i = 1,\ldots, N-1,\\[7pt] f_N = \dfrac{\mu_2 + \eta_2\beta_N}{1 - \alpha_N\eta_2},\\ f_i = \alpha_{i+1}f_{i+1} + \beta_{i+1}, i = N-1,\ldots, 0. \end{cases} $$

Моя программа считывает значения $N$, $M$, $h$, $\tau$, $\eta_1$, $\mu_1$, $\eta_2$, $\mu_2$, $u$, момент $t_k$, в который будет строиться график и функцию $\varphi(x, 0)$ (начальное распределение).

Ниже --- код проги (fi[N+1][M+1] --- сетка; fi_x_0(x) --- $\varphi(x, 0)$). Идея: заполняется нижний ряд точек (fi[i][0]), а для каждого следующего ряда применяем метод прогонки, беря значения из предыдущего.

  def solve
    # filling fi(x, 0)
    for i in 0 .. @n
      @fi[i][0] = fi_x_0( i*@h )
    end
    # tdma
    alpha, beta = [], []
    alpha[1], beta[1] = @eta_1, @mu_1
    aa_i = -@tau * @u / (2.0 * @h)
    bb_i =  @tau * @u / (2.0 * @h)
    cc_i = -1.0
    for j in 0 ... @m
      for i in 1 ... @n
        alpha[i+1] = bb_i / (cc_i - alpha[i] * aa_i)
        beta[i+1]  = (aa_i * beta[i] + @fi[i][j]) / (cc_i - alpha[i] * aa_i)
      end
      @fi[@n][j+1] = (@mu_2 + @eta_2 * beta[@n]) / (1 - alpha[@n] * @eta_2)
      (@n - 1).downto( 0 ) do |i|
        @fi[i][j+1] = alpha[i+1] * @fi[i+1][j+1] + beta[i+1]
      end
    end
  end

Но на $N = 100$, $M = 100$, $h = 0.01$, $\tau = 0.01$, $\eta_1 = 0.5$, $\mu_1 = 0.2$, $\eta_2 = 0.3$, $\mu_2 = 0.6$, $u = 1$, $t_k = 0.01$,

$$ \varphi(x, 0) = \begin{cases} 0.5\sin(8\pi x - \pi/2) + 1.5, & x \leqslant 0.25,\\ 1, & \text{ else } \end{cases} $$

получается такая фигня:

а хочется получить что-то такое (пардон за качество):


У меня какая-то ошибка в идее реализации алгоритма или в начальных данных. Буду признателен за любую помощь)

[Imperator]

Ошибка может сидеть где угодно. Проверьте для начала отдельно на какой-то системе с трехдиагональной матрицей работоспособность метода прогонки.
Затем проверьте, правильно ли собирается СЛАУ (матрица системы), которую требуется решить.
Рассмотрите для начала какую-нибудь небольшую сетку (например, когда $i=0,...,3$, $j=0,...,3$ (16 узлов) или еще меньше).

[worm2]

Что-то схема, да и сама постановка задачи, вызывают подозрения.
Я, конечно, не великий специалист, но тот факт, что Вам приходится задавать условие на правой границе:

Цитата:

Граничные условия задаём в виде ... $f_N = \eta_2f_{N-1} + \mu_2$

настораживает.
Для уравнения переноса граничное условие задаётся только на одной границе: правой, если скорость $u > 0$, и левой, если меньше нуля.
Причём это условие должно быть уже в исходной постановке, до того, как мы начинаем строить разностные схемы для его решения (хотя бы потому, что задачу неплохо бы корректно поставить, прежде чем пытаться её решать).
Т.е. кроме начального условия $\varphi(x,0)=\dots$ нужно ещё условие $\varphi(0,t)=\dots$ (здесь предполагаю, что $u > 0$), однако я такового не вижу.
Но даже после того, как задача будет корректно поставлена, схема, требующая задания дополнительного (отсутствующего в исходной постановке) граничного условия, априори вызывает сомнения. Я, конечно, понимаю, что хочется получить второй порядок аппроксимации по $x$, но не вижу, как эта идея могла бы здесь прокатить. Кто Вам указал использовать именно такую схему?

[Pete[r]]

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

Задание: используя неявную разностную схему $\dfrac{\varphi_i^{j+1} - \varphi_i^j}{\tau} + u\dfrac{\varphi_{i+1}^{j+1} - \varphi_{i-1}^{j+1}}{2h} = 0$, найти числ. реш-е ур-я конвект. переноса в классе сеточных функций $\varphi_i^j$, определённом на (нарисованном выше) 4-хклеточном шаблоне на отрезке $0\leqslant x = ih\leqslant s = Nh$ в промежутке времени $0\leqslant t = j\tau\leqslant M\tau$, если начальное распределение искомой функции описывается заданным преподавателем соотношением. Для получения решения использовать метод прогонки. Решение получить для значений сеточного параметра $\gamma < 1$, $\gamma = 1$, $\gamma > 1$ [надо проверить, что схема даёт устойчивые решения при любых $\gamma$].

Дело в том, что я численными методами занимался довольно давно и совсем немного, так что мне самому непонятен физический смысл констант $\eta_i$ и $\mu_i$. Т.е. я понимаю, как на их основе решается трёхдиагональная система, но в чём проблема в данном случае, не осознаю.
Ниже, на средней картинке, написано, что $\eta_i$ определяют 1-й, 2-ой или 3-й род граничных условий, но мне это мало о чём говорит...

[worm2]

Чё-та какая-то попсовая методичка...

В формуле (2.23) ошибка: написано $F_i=\phi_i^n$, а должно быть $F_i=-\phi_i^n$.

[Pete[r]]

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

Чё-та какая-то попсовая методичка...

Что ж, как и многие другие.

Даже если взять $F_i=-\phi_i^n$, график просто перевернётся. Ладно, спасибо за советы, буду проверять все формулы...

[worm2]

Перевернётся и станет почти правильным.
Дело в том, что в этой методичке ещё и в рисунке ошибка.
В начальном значении $\varphi(x,0)$ максимум будет в точке 0.125, соответственно, в решении он чуть сдвинется вправо.
А в отсканированном рисунке этот максимум в районе 0.25, что неправильно.
А то, что у Вас левый хвост уходит ниже значения 1, это вполне объяснимо. Ведь у Вас $f_0 = \eta_1f_1 + \mu_1$, что является весьма странным условием, коэффициенты просто взяты от балды (так же, как и коэффициенты для правого краевого условия).
Почему взяты от балды, конечно, понятно. Неоткуда взять правильное граничное условие при $x=0$, ибо у исходной задачи оно не задано. Тут всё что угодно можно напридумывать.

Я, впрочем, подозреваю, что подразумевалось что-то вроде того, что начальное условие задано на всём $\mathbb R$, например, так:

$$ \varphi(x, 0) = \begin{cases} 0.5\sin(8\pi x - \pi/2) + 1.5, & 0 \leqslant x \leqslant 0.25,\\ 1, & \text{ else } \end{cases} $$

(я добавил $0 \leqslant x$), ну и, типа, решается не краевая задача, а задача Коши.
Тогда, действительно, граничное условие легко восстанавливается:
$\varphi(0, t)=0$, откуда и условие на $f_0$ тривиально выписывается, и никаких $\eta_1$, $\mu_1$ не нужно. Но от $\eta_2$, $\mu_2$ избавиться с такой схемой, вроде, не получится, ввиду чего, при задании их от балды, и получаются "глюки" на правом конце отрезка.
Так что меня в Вашем решении ничего не удивляет, весьма вероятно, что Вы запрограммировали всё правильно.

Добавлено спустя 11 минут 19 секунд:

А, собственно, с чего я взял, что решается краевая задача, а не задача Коши? Видимо, задача Коши и решается
Меня сбило с толку условие $f_0 = \eta_1f_1 + \mu_1$ для уже разностной задачи.

Всё равно, в методичке как-то очень неаккуратно... Похоже, $\varphi(x,0)$ задано именно так, как я написал, сужу по отсканированному рисунку в методичке, иначе решение при уменьшении x меньше нуля должно было расти, а оно остаётся равым 1.

Всё вышенаписанное остаётся в силе, за исключением замечания, что нужно задавать $\varphi(0,t)$. Не нужно. Тем не менее, формулы все остаются (в т.ч. и правильное условие на левом конце сетки: $f_0=0$).

[Pete[r]]

Спасибо большое за подробное объяснение!

Надеюсь, так и подразумевалось, впрочем, от меня требуется только реализовать алгоритм в виде программы, а учитывая лажовость методички, никакого желания тратить драгоценное время на поиск опечаток нет