Разностная схема для уравнения переноса. Монотонность

[Challenger]

Здравствуйте! Помогите пожалуйста с задачей:
Дано уравнение переноса:
$u_t+cu_x = f, ~ c = const > 0$
Нужно составить разностную схему с порядком аппроксимации $O(\tau, h^2)$ на 4-х точечном шаблоне - квадрат (не стал вставлять картинку, думаю понятно). Исследовать полученную схему на устойчивость и монотонность (по Фридрихсу).
Значит, схема будет следующей (метод прямых):
$\frac{1}{2}\Bigl[\frac{y_m^{n+1} - y_m^n}{\tau} + \frac{y_{m+1}^{n+1} - y_{m+1}^n}{\tau} \Bigr] + \frac{c}{2}\Bigl[\frac{y_{m+1}^{n+1} - y_m^{n+1}}{h} +\frac{y_{m+1}^n - y_m^n}{h} \Bigr] = \varphi_m$
$m$ - нумерует узлы по $x$, $n$ - по $t$.
Устойчивость я вроде как показал, получилось, что она безусловно устойчива.
Но вот с монотонностью по Фридрихсу проблемы, так как схема неявная, а я не нашел достаточного условия монотонности для неявных схем. Единственное, что Калиткин "Численные методы, 1978" на стр. 348 пишет, что
"если двухслойная однородная схема неявна, то ее можно преобразовать к явной форме $\hat{y}_m = \sum\limits_{l} \beta_l y_{m+l}$, где пределы суммы по $l$ бесконечны, и затем применить признак монотонности.
Но я не в силах понять, что значит "пределы суммы по $l$ бесконечны", ведь во-первых, у нас фиксирован шаблон, во-вторых чем это мне поможет? То есть какие коэффициенты будут стоять перед другими $y_{m+l}$?

[vasya321]

Здравствуйте. А как мы будем выполнять расчет по данной схеме? У нас известны какие-то граничные условия на левой границе по $x$: $u(0,t) = g(t)$, а также начальное распределение $u$: $u(x,0) = h(x)$. То есть, нам известны $u^{n}_{0} = g(t_n), n=0,1, \ldots, N_{t}$ и $u^{0}_{m} = h(x_m), m=0,1, \ldots, N_{x}$, здесь $N_{t}, N_{x}$ — число узлов сетки по $t$ и $x$.

Таким образом, в процессе расчета нам всегда известны $u^{n}_{m}$, $u^{n+1}_{m}$, $u^{n}_{m+1}$ по которым рассчитывается $u^{n+1}_{m+1}$

$$ \begin{tikzpicture} \draw (0,0) -- (1,0) -- (1,1) -- (0,1) -- (0,0); \filldraw (0,0) circle (0.1) node[left] {$u^{n}_{m}$}; \filldraw (1,0) circle (0.1) node[right] {$u^{n}_{m+1}$}; \filldraw[fill=white] (1,1) circle (0.1) node[right] {$u^{n+1}_{m+1}$}; \filldraw (0,1) circle (0.1) node[left] {$u^{n+1}_{m}$}; \end{tikzpicture} $$

Для доказательства монотонности необходимо выразить значение $u^{n+1}_{m+1}$ через известные. Если хотя бы один из коэффициентов при $u^{n}_{m}$, $u^{n+1}_{m}$, $u^{n}_{m+1}$ отрицательный, то схема немонотонна по теореме Годунова о сохранении монотонности (Theorem 1: Monotonicity preserving)

Кстати, данная схема имеет такой же порядок точности по времени как и по координате.

[Challenger]

Так я совсем запутался... Выходит, это схема явная, не смотря на то, что на верхнем слое два узла?

vasya321 в сообщении #1016538 писал(а):

Таким образом, в процессе расчета нам всегда известны $u^{n}_{m}$, $u^{n+1}_{m}$, $u^{n}_{m+1}$ по которым рассчитывается $u^{n+1}_{m+1}$

А мы не считаем $u^{n+1}_{m}$ неизвестным?

vasya321 в сообщении #1016538 писал(а):

Кстати, данная схема имеет такой же порядок точности по времени как и по координате.

Странно, у меня в конспекте написано $O(\tau, h^2)$.

[vasya321]

Challenger в сообщении #1016571 писал(а):

Выходит, это схема явная, не смотря на то, что на верхнем слое два узла?

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

Думаю, что Калиткин под "пределы суммы по $l$ бесконечны" имеет в виду, что для неявных линейных схем каждое значение на неизвестном слое зависит от всех значений на известном. Например, рассмотрим неявную схему уравнения диффузии

$$ \frac{u^{n+1}_m - u^{n}_m}{\Delta t} = \frac{u^{n+1}_{m-1} - 2 u^{n+1}_{m} + u^{n+1}_{m-1}}{\Delta x^2} $$

Как известно, эта схема приводит к системе трех точечных уравнений вида

$$ \beta_{m-1,m} u^{n+1}_{m-1} + \beta_{m,m} u^{n+1}_{m} + \beta_{m-1,m} u^{n+1}_{m-1} = f_{m}(u^{n}_m), \quad m=1,2 \ldots, N_x-1, $$

которую можно переписать в матричном виде

$$ B U^{(n+1)} = F(U^{(n)}), $$

где $B$ - матрица коэффициентов $\beta$, $U^{(n+1)}$ - искомый вектор решения на неизвестном временном слое $u^{n+1}_m$.

Решение $U^{(n+1)}$ получаем как

$$ U^{(n+1)} = B^{-1} F(U^{(n)}) = C U^{(n)} $$

Если в матрице $B$ много нулей, то матрица $B^{-1}$ полностью заполнена. Поскольку $u^{n}_m$ входит в $f_m$ линейно, то коэффициенты перед $u^{n}_m$ можно вынести в отдельную матрицу и, перемножив с $B^{-1}$, получить матрицу $C$. Мы теперь получили как бы явную схему, поскольку:

$$ u^{n+1}_{m} = \sum \limits_{l=0}^{N_x} \gamma_{l,m} u^{n}_l, $$

где $\gamma_{l,m}$ - это элементы матрицы $C$. Условие монотонности теперь выражается в положительности всех $\gamma_{l,m}$.

Challenger в сообщении #1016571 писал(а):

Странно, у меня в конспекте написано $O(\tau, h^2)$.

Вы можете показать порядок точности, разложив $u(x,t)$ в ряд Тейлора для двух переменных в окрестности центра квадрата, и найти выражения для производных с помощью этого ряда. Получается, что производные в уравнении $u_t + c u_x=0$ аппроксимируются центральными конечными разностями в центре квадрата.

[Challenger]

vasya321 в сообщении #1018954 писал(а):

Думаю, что Калиткин под "пределы суммы по $l$ бесконечны" имеет в виду, что для неявных линейных схем каждое значение на неизвестном слое зависит от всех значений на известном. Например, рассмотрим неявную схему уравнения диффузии...

Насколько я понял, вы просто описали метод решения неявных разностных схем. А обращать матрицу не сложно будет?
Оказывается, Калиткин там же (на стр. 349) приводит пример как сводить неявную схему к явной. Но там все просто и понятно, в отличие от моего случая.

Challenger в сообщении #1016571 писал(а):

Странно, у меня в конспекте написано $O(\tau, h^2)$.

vasya321 в сообщении #1018954 писал(а):

Вы можете показать порядок точности, разложив $u(x,t)$ в ряд Тейлора для двух переменных в окрестности центра квадрата, и найти выражения для производных с помощью этого ряда. Получается, что производные в уравнении $u_t + c u_x=0$ аппроксимируются центральными конечными разностями в центре квадрата.

Да, пардон, не туда посмотрел. Вы абсолютно правы. Тогда эта схема нам не подходит, нужно обязательно $O(\tau, h^2)$. Думаю надо просто изменить разностную производную по $x$

$$\frac{1}{2}\Bigl[\frac{y_m^{n+1} - y_m^n}{\tau} + \frac{y_{m+1}^{n+1} - y_{m+1}^n}{\tau} \Bigr] + c\Bigl[\frac{y_{m+1}^{n+1} - y_m^{n+1}}{h} \Bigr] = \varphi_m ?$$

[vasya321]

Challenger в сообщении #1021047 писал(а):

нужно обязательно $O(\tau, h^2)$. Думаю надо просто изменить разностную производную по $x$

$$\frac{1}{2}\Bigl[\frac{y_m^{n+1} - y_m^n}{\tau} + \frac{y_{m+1}^{n+1} - y_{m+1}^n}{\tau} \Bigr] + c\Bigl[\frac{y_{m+1}^{n+1} - y_m^{n+1}}{h} \Bigr] = \varphi_m ?$$

Да, всё верно.

Challenger в сообщении #1021047 писал(а):

Оказывается, Калиткин там же (на стр. 349) приводит пример как сводить неявную схему к явной. Но там все просто и понятно, в отличие от моего случая.

Почему вы не можете использовать такой же приём как Калиткин? Выразить $y^{n+1}_{m+1}$ через другие значения $y$ расчетной схемы, подставить в эту формулу $y^{n+1}_{m}$, проследить закономерность коэффициентов при $y^{n}_{m+1-l}, \; l=0,1 \ldots N_x-1$ и получить формулу для явной схемы. Надо помнить, что $y^{n+1}_{0}$ известно из граничных условий.

[Challenger]

vasya321 в сообщении #1021274 писал(а):

Почему вы не можете использовать такой же приём как Калиткин? Выразить $y^{n+1}_{m+1}$ через другие значения $y$ расчетной схемы, подставить в эту формулу $y^{n+1}_{m}$, проследить закономерность коэффициентов при $y^{n}_{m+1-l}, \; l=0,1 \ldots N_x-1$ и получить формулу для явной схемы. Надо помнить, что $y^{n+1}_{0}$ известно из граничных условий.

Да, с этим наконец-то справился!

vasya321 в сообщении #1021274 писал(а):

Challenger в сообщении #1021047 писал(а):

нужно обязательно $O(\tau, h^2)$. Думаю надо просто изменить разностную производную по $x$

$$\frac{1}{2}\Bigl[\frac{u_m^{n+1} - u_m^n}{\tau} + \frac{u_{m+1}^{n+1} - u_{m+1}^n}{\tau} \Bigr] + c\Bigl[\frac{u_{m+1}^{n+1} - u_m^{n+1}}{h} \Bigr] = \varphi_m ?$$

Да, всё верно.

Хорошо, но теперь меня смущает порядок аппроксимации: если раскладывать по формуле Тейлора в $(t_n +\frac{\tau}{2}, x_m + \frac{h}{2})$, то по $\tau$ получается 2-ой порядок. Смотрите
$u_{m+1}^{n+1} = u_{m+\frac{1}{2}}^{n+\frac{1}{2}}+ \frac{\tau}{2}u_t+\frac{\tau^2}{8}u_{tt}+O(\tau^3),$
$u_{m+1}^{n} = u_{m+\frac{1}{2}}^{n+\frac{1}{2}}- \frac{\tau}{2}u_t+\frac{\tau^2}{8}u_{tt}+O(\tau^3)$
Соответственно,
$\frac{u_{m+1}^{n+1} - u_{m+1}^n}{\tau} = u_t + O(\tau^2)$
Аналогично для $\frac{u_{m}^{n+1} - u_{m}^n}{\tau}$
Где я ошибаюсь?
Плюс еще один глупый вопрос: почему вообще можно раскладывать по ф. Тейлора в точке, которой нет в шаблоне?
И еще: что делать с аппроксимацией правой части $f$? Не понизит ли она случаем порядок $O(\tau, h^2)$?

[vasya321]

Challenger в сообщении #1021451 писал(а):

Хорошо, но теперь меня смущает порядок аппроксимации: если раскладывать по формуле Тейлора в $(t_n +\frac{\tau}{2}, x_m + \frac{h}{2})$, то по $\tau$ получается 2-ой порядок. Смотрите
$u_{m+1}^{n+1} = u_{m+\frac{1}{2}}^{n+\frac{1}{2}}+ \frac{\tau}{2}u_t+\frac{\tau^2}{8}u_{tt}+O(\tau^3),$
$u_{m+1}^{n} = u_{m+\frac{1}{2}}^{n+\frac{1}{2}}- \frac{\tau}{2}u_t+\frac{\tau^2}{8}u_{tt}+O(\tau^3)$
Соответственно,
$\frac{u_{m+1}^{n+1} - u_{m+1}^n}{\tau} = u_t + O(\tau^2)$
Аналогично для $\frac{u_{m}^{n+1} - u_{m}^n}{\tau}$
Где я ошибаюсь?

Похоже, вы забываете, что на значение $u^{n+1}_{m+1}$ влияют также производные по $x$. Расположим систему координат в центре квадрата и разложим $u(x,t)$ в ряд Тейлора для двух переменных в окрестности (0,0) [раскладывать надо до 3-х частных производных]:

$$ u(x,t) = u(0,0) + x \frac{\partial u}{\partial x} + t \frac{\partial u}{\partial t} + \frac{1}{2 !} \left[ x^2 \frac{\partial^2 u}{\partial x^2} + 2 x t \frac{\partial u}{\partial x \partial t} + t^2 \frac{\partial^2 u}{\partial t^2} \right] + \cdots $$

Теперь подставим значения $u^{n+1}_{m+1} = u(\Delta x /2,\Delta t /2)$, $u^{n}_{m+1} = u(\Delta x /2, -\Delta t /2)$ и так далее в к-р аппроксимацию.

Challenger в сообщении #1021451 писал(а):

почему вообще можно раскладывать по ф. Тейлора в точке, которой нет в шаблоне?

Какие есть конкретные препятствия этому?

Challenger в сообщении #1021451 писал(а):

И еще: что делать с аппроксимацией правой части $f$?

Чтобы уравнение $u_t + c u_x = f$ точно выполнялось в точке 0,0 надо взять $f$ в этой точке, т.е. $\varphi_{m+\frac12}^{n+\frac12}$

Challenger в сообщении #1021451 писал(а):

Не понизит ли она случаем порядок $O(\tau, h^2)$?

Подставьте в формулу для $f$ разложения в ряд Тейлора $u(x,t)$, а саму $f$ в к-р аппроксимацию и посмотрите.

[Challenger]

vasya321 в сообщении #1021547 писал(а):

Похоже, вы забываете, что на значение $u^{n+1}_{m+1}$ влияют также производные по $x$. Расположим систему координат в центре квадрата и разложим $u(x,t)$ в ряд Тейлора для двух переменных в окрестности (0,0) [раскладывать надо до 3-х частных производных]:

$$ u(x,t) = u(0,0) + x \frac{\partial u}{\partial x} + t \frac{\partial u}{\partial t} + \frac{1}{2 !} \left[ x^2 \frac{\partial^2 u}{\partial x^2} + 2 x t \frac{\partial u}{\partial x \partial t} + t^2 \frac{\partial^2 u}{\partial t^2} \right] + \cdots $$

Теперь подставим значения $u^{n+1}_{m+1} = u(\Delta x /2,\Delta t /2)$, $u^{n}_{m+1} = u(\Delta x /2, -\Delta t /2)$ и так далее в к-р аппроксимацию.

Ага, вас понял. Теперь действительно получается. В производной по $x$ остается $\tau$ из-за смешанной производной $u_{xt}$.

Challenger в сообщении #1021451 писал(а):

почему вообще можно раскладывать по ф. Тейлора в точке, которой нет в шаблоне?

vasya321 в сообщении #1021547 писал(а):

Какие есть конкретные препятствия этому?

Да, наверное никаких, глупость сморозил.Спасибо!