Начально-краевая задача для уравнения переноса.

[a3ro]

Нужно решить численно начально-краевую задачу для уравнения переноса $u_t+u_x=0$
Начальное условие: $u(0,x) = 1, x\leqslant5 $ и $u(0,x) = 0, 5\leqslant x \leqslant20$ и граничное условие $u(t,0) = 1$.
Нужно использовать начальную сетку из 41 точки при $h = 0,5$ и провести расчёт до значения $T = 9$. Решить для чисел Куранта $\tau/h = 1.0, 0.6, 0.3$. Рекомендованная разностная схема Лакса-Вендроффа.

Моё решение:
Схема Лакса-Вендроффа: $\frac{u_j^{i+1} - u_j^i}{\tau} + a \frac{u_{j+1}^i - u_{j-1}^i}{2h} - \frac{a^2\tau}{2}\frac{u_{j+1}^i - 2u_j^i + u_{j-1}^i}{h^2} = 0 $
Выражаю $u_j^{i+1}$:
$u_j^{i+1} = u_j^i + \frac{a^2\tau^2}{2h^2}(u_{j+1}^i - 2u_j^i + u_{j-1}^i) - \frac{a\tau}{2h}(u_{j+1}^i - u_{j-1}^i)$
В нашем случае (для начала число куранта будет 1.0) $a = 1$, $\tau = 0.5$:
$u_j^{i+1} = u_j^i + \frac{1}{2}(u_{j+1}^i - 2u_j^i + u_{j-1}^i) - \frac{1}{2}(u_{j+1}^i - u_{j-1}^i)$
Ну и тут получается бред, что $u_j^{i+1} = u_{j-1}^i$ (*)

В программе для начала, чтобы не путаться, вбил статический массив (с учётом шага по вертикали и горизонтали по 0.5):

double T = 9;
   double L = 20;
   double u[19][41];

Дальше задавал граничные и начальные условия:

   for (int i = 0; i <= 18; i++){
      u[i][0] = 1;
   }
   for (int j = 0; j <= 10; j++){
      u[0][j] = 1;
   }
   for (int j = 11; j <= 40; j++){
       u[0][j] = 0;
   }

Ну и для следующего момента времени функцию u выражал по вышеприведенной формуле:

for (int i = 0; i <= 17; i++){
      for (int j = 1; j <= 39; j++) {
         u[i + 1][j] = u[i][j] + 1 / 2 * (u[i][j + 1] - 2 * u[i][j] + u[i][j - 1]) - 1 / 2 * (u[i][j + 1] - u[i][j - 1]);
      }
   }

И вот, собственно, результат. Очевидно, это не то.
При числах куранта 0.6 и 0.3 получается тоже самое, хотя формула (*) уже не выглядит такой тривиальной.


Я понимаю, что ошибся, и что ошибка моя в самом начале и связана, как мне кажется, с недопониманием. Прошу навести на верный путь .

[worm2]

a3ro писал(а):

Ну и тут получается бред, что $u_j^{i+1} = u_{j-1}^i$ (*)

Почему бред? Так и должно было получиться.

a3ro писал(а):

Используется синтаксис C

for (int i = 0; i <= 17; i++){
  for (int j = 1; j <= 39; j++) {
    u[i + 1][j] = u[ i ][j] + 1 / 2 * (u[ i ][j + 1] - 2 * u[ i ][j] + u[ i ][j - 1]) - 1 / 2 * (u[ i ][j + 1] - u[ i ][j - 1]);
  }
}

Не стесняйтесь, пишите, как получилось после упрощений:

Используется синтаксис C

    u[i + 1][j] = u[ i ][j - 1];

Здесь всё верно. Вот результат уже немного смущает. Единички должны на каждом шаге на одну вверх смещаться. Точнее, вдоль прямых $x-t=\rm{const}$ решение должно быть постоянным: где на пересечении этой прямой с границей были нулики, там и будет 0, а где были единички, будет 1. Это верно для точного решения задачи. Но в данном случае (когда число Куранта = 1) это верно и для сеточного решения, ибо здесь сеточное совпадает с точным (в тех точках, в которых вычисляется).
Что-нибудь понятно из моего объяснения? Ключевой вопрос, который нужно здесь понять — вопрос характеристик уравнения переноса. Знаете, что это такое?

[a3ro]

Цитата:

Единички должны на каждом шаге на одну вверх смещаться. Точнее, вдоль прямых $x-t=\rm{const}$ решение должно быть постоянным

Если учитывать, что у меня ось $\tau$ направлена "вниз", то наверное в каждой следующей строке должно быть на одну единичку больше. Т.е. от "единиц" идут "единицы" вправо вниз, от нулей должны идти нули. Поправьте, если сказал не так.
Значит ошибка всё-таки была не в начале, а в коде, сейчас попробую пофиксить.

Цитата:

Ключевой вопрос, который нужно здесь понять — вопрос характеристик уравнения переноса. Знаете, что это такое?

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

[worm2]

Ну да, всё правильно. Кроме ответа, почему-то.
Возможно, компилятор воспринимает 1 / 2 как целочисленное деление и выдаёт в промежуточном результате 0. Это ж Си, он такое любит. Я бы во всех таких местах поставил 0.5. Или 1.0 / 2.0.

[a3ro]

Пофиксил. Я уже забыл о куче подводных камней в языке Си да и в программировании в целом.


Получается для числа Куранта $q=1$ сгущение сеток для данных граничных и начальных условий не будет влиять на погрешность численного решения (т.к. погрешности нет). Но для чисел Куранта $q_2=0.6$ и $q_3=0.3$ такая погрешность появится.
И дальше немного глупый вопрос, который вытекает из условия задачи: "просят оценить порядок сходимости схемы сравнением изменения погрешности численного решения по сравнению с точным на последовательности сгущающихся сеток при фиксированном числе Куранта (применить сгущение сеток вдвое)".
Т.е. порядок такой:
1) Мы берем число Куранта. (А дальше у меня в голове путаница...)
2) Для каждого параметра h (0.5, 0.25 или 0.125) сетки мы берем разность $([u]_h - u)$, где u точное решение, и сравниваем его c $h^k$ $(ch^k ?)$ где k и есть наш порядок сходимости.
Т.е. $([u]_h - u) \leqslant h^kc$.
У нас такое соотношение должно выполняться для каждой точки? Как находить число c?

[worm2]

a3ro писал(а):

Как находить число c?

Тут главное оценить порядок (как вас и просят в условии), т.е. $k$. Например, находите максимум $([u]_h - u)$ по всем точкам, а потом то же самое для уменьшенного в 2 (3, 4, 10, ...) раз шага. Максимум погрешности должен при этом уменьшиться примерно в $2^k$ ($3^k$, $4^k$, $10^k$, ...) раз. После этого и $c$ можно оценить (хотя этого и не просят).

[a3ro]

Собственно решил для примера вывести решение для сетки $h = 0.5, q = 0.6$

(Пришлось в 2 картинки из-за ограничения по ширине).
Где приведенная численная схема выглядит: $u_j^{i+1} = 0.64u_j^i + 0.48u_{j-1}^i - 0.12u_{j+1}^i$

Ну и вот, что я подумал. Тут у нас полный диапазон значений от 0 до 1.2. А точные решения в данной задаче либо 1, либо 0. Максимум $([u]_h-u)$ по всем точкам будет 0.5 (как видно из скрина). Я сделал тоже самое для вдвое сгущенной сетки ($h = 0.25$), там тоже встречаются значения $u_h = 0.5$, а значит и максимум $([u]_h-u)$ будет тоже 0.5.

Поправьте, пожалуйста, если я ошибаюсь.

[worm2]

a3ro писал(а):

$u_j^{i+1} = 0.64u_j^i + 0.48u_{j-1}^i - 0.12u_{j+1}^i$

Ниже неправильные рассуждения:
Глядя только на это, без предыдущих выкладок, можно сразу сказать, что условие Куранта тут нарушено. Число Куранта > 1. В результате возник отрицательный коэффициент (-0.12) там, где (при выполнении условия Куранта) должны были быть только положительные.
Это приводит к нарушению устойчивости схемы и делает её бесполезной чуть менее, чем полностью. По результату это тоже видно, как только значения начинают превышать максимум краевых значений 1.0: при устойчивой схеме такое немыслимо Думаю, дальше там и отрицательные значения повылезают. Ошибки округления растут экспоненциально, и довольно быстро забивают всё решение. Нет смысла говорить о порядке точности в таких условиях.

[a3ro]

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

a3ro писал(а):

$u_j^{i+1} = 0.64u_j^i + 0.48u_{j-1}^i - 0.12u_{j+1}^i$

Глядя только на это, без предыдущих выкладок, можно сразу сказать, что условие Куранта тут нарушено. Число Куранта > 1. В результате возник отрицательный коэффициент (-0.12) там, где (при выполнении условия Куранта) должны были быть только положительные.
Это приводит к нарушению устойчивости схемы и делает её бесполезной чуть менее, чем полностью. По результату это тоже видно, как только значения начинают превышать максимум краевых значений 1.0: при устойчивой схеме такое немыслимо Думаю, дальше там и отрицательные значения повылезают. Ошибки округления растут экспоненциально, и довольно быстро забивают всё решение. Нет смысла говорить о порядке точности в таких условиях.

Так, давайте распишу тогда подробнее:

Цитата:

Вот наша формула:
$u_j^{i+1} = u_j^i + \frac{a^2\tau^2}{2h^2}(u_{j+1}^i - 2u_j^i + u_{j-1}^i) - \frac{a\tau}{2h}(u_{j+1}^i - u_{j-1}^i)$

$a = 1, q=\tau/h = 0.6$. Тогда $\tau^2/h^2 = 0.36. Подставляю в уравнение выше:
$u_j^{i+1} = u_j^i + 1/2 \cdot 0.36 \cdot (u_{j+1}^i - 2u_j^i + u_{j-1}^i) - 1/2 \cdot 0.6 \cdot(u_{j+1}^i - u_{j-1}^i)$
$u_j^{i+1} = u_j^i + 0.18 \cdot (u_{j+1}^i - 2u_j^i + u_{j-1}^i) - 0.3 \cdot (u_{j+1}^i - u_{j-1}^i)$
И вот тут у нас вылазит минус перед членом $u_{j+1}^i$. Для числа Куранта 0.3 по таким же рассуждениям тоже появится отрицательный член.

[worm2]

Да, я наврал в предыдущем сообщении, поспешил Хотел спинным мозгом обойтись, а тут надобно головным думать

Условие Куранта выполняется. И формулы у вас правильные. И отрицательные коэффициенты для схемы Лакса–Вендроффа неизбежны. Но устойчивость всё-таки есть, просто она более хитро проверяется, не на основе принципа максимума, а по фон Нейману.

Соответственно, и норму невязки придётся считать не как максимум $([u]_h - u)$ по всем точкам, а как-то усредняя погрешность по всем точкам.
Например, в сеточном аналоге нормы $L_2$ как-то так:

$$||u_h - u||=\sqrt{\frac{1}{(N_x+1) (N_t+1)}\sum\limits_{i=0}^{N_t}\sum\limits_{j=0}^{N_x}\left(u_j^i-u_{\text{точн}}(x_j, t_i)\right)^2}.$$

И вот эта погрешность должна убывать как $k$-я степень шага при уменьшении шагов с сохранением условия Куранта.

[a3ro]

Благодарю, вчера еще увидел, но уснул на 15 часов

Не будет ли в данном случае халтурой, если я буду просто каждый элемент численных решений сравнивать с 0 и 1 и брать наименьшую разность, чтобы подставить в фон Неймана?

Ну и вопрос (возможно не совсем по теме) по Лаксу-Вендроффу: как посроить эту схему для стационарного уравнения Эйлера:
$v\frac{dv}{dx}$ = -$\frac{1}{\rho}$$\frac{dp}{dx}$.
Моё непонимание в том, что тут 2 функции, а дифференцируются по 1 переменной, а не наоборот.
Или в данных случаях переворачивают равенство (хотя можно ли это делать, когда у нас скорость и давление функции нескольких переменных для 2D случая), накладывая условия для области определения?

[worm2]

a3ro в сообщении #1315917 писал(а):

Не будет ли в данном случае халтурой, если я буду просто каждый элемент численных решений сравнивать с 0 и 1 и брать наименьшую разность, чтобы подставить в фон Неймана?

Не знаю. Боюсь, не получится.

a3ro в сообщении #1315917 писал(а):

Моё непонимание в том, что тут 2 функции, а дифференцируются по 1 переменной, а не наоборот.

Когда одна переменная, а функций две, это гораздо проще. Пара обычных дифуров, не в частных производных. Можно обойтись без разностных схем, а методами вроде Рунге-Кутта. Только проблема, что неизвестных больше, чем уравнений Ну то есть одного уравнения Эйлера мало. Или вам давление дано, и вы по нему можете вычислить скорость, или наоборот. Или ещё что-то.
Когда переменных две или больше, тут уже будут уравнения в частных производных, и схема понадобится. Только та же проблема остаётся: не хватает одного этого уравнения, нужно ещё что-то, либо давление, либо скорость, либо ещё какое-то уравнение. Нужно поработать над постановкой, а потом уже задумываться над методами решения.
И потом, тут вы перескакиваете от модельной математической задачи к более сложной задаче гидродинамики, не слишком ли быстро?

[a3ro]

Цитата:

Не знаю. Боюсь, не получится.

Да просто не хочу "вручную" забивать точные данные для каждой точки в сетках, где число куранта не 0, чтобы сравнивать (хочу обойтись спинным мозгом )

Цитата:

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

Я опечатался в предыдущем сообщении, там p - это импульс.
Задача заключается в том, что нужно в рамках краевой задачи для уравнения Эйлера решить задачу об обтекании кругового цилиндра поперек потока. Ничего больше не дали
Тут нужно еще уравнение $divV = 0$ - уравнение "неразрывности" для жидкости (плотность $\rho$ const).

Цитата:

И потом, тут вы перескакиваете от модельной математической задачи к более сложной задаче гидродинамики, не слишком ли быстро?

Ну гидродинамику я немного знаю, и подобную задачу решал через комплексный потенциал ,находил распределение скорости на поверхности (закон от угла $V = 2V_\infty \sin\theta$, где $\theta$ - угол "от центра"). Но с уравнением Эйлера дел на конкретных примерах не имел.

[worm2]

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

Кстати, вроде бы $p$ — всё-таки давление.

[worm2]

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