2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу 1, 2  След.
 
 Начально-краевая задача для уравнения переноса.
Сообщение26.05.2018, 13:40 
Нужно решить численно начально-краевую задачу для уравнения переноса $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 получается тоже самое, хотя формула (*) уже не выглядит такой тривиальной.
Изображение

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

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение26.05.2018, 19:37 
Аватара пользователя
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) это верно и для сеточного решения, ибо здесь сеточное совпадает с точным (в тех точках, в которых вычисляется).
Что-нибудь понятно из моего объяснения? Ключевой вопрос, который нужно здесь понять — вопрос характеристик уравнения переноса. Знаете, что это такое?

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

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

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

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

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

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение26.05.2018, 22:34 
Пофиксил. Я уже забыл о куче подводных камней в языке Си да и в программировании в целом.
Изображение

Получается для числа Куранта $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?

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение26.05.2018, 22:49 
Аватара пользователя
a3ro писал(а):
Как находить число c?
Тут главное оценить порядок (как вас и просят в условии), т.е. $k$. Например, находите максимум $([u]_h - u)$ по всем точкам, а потом то же самое для уменьшенного в 2 (3, 4, 10, ...) раз шага. Максимум погрешности должен при этом уменьшиться примерно в $2^k$ ($3^k$, $4^k$, $10^k$, ...) раз. После этого и $c$ можно оценить (хотя этого и не просят).

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение27.05.2018, 17:45 
Собственно решил для примера вывести решение для сетки $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.

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

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение28.05.2018, 00:21 
Аватара пользователя
a3ro писал(а):
$u_j^{i+1} = 0.64u_j^i + 0.48u_{j-1}^i - 0.12u_{j+1}^i$
Ниже неправильные рассуждения:
Глядя только на это, без предыдущих выкладок, можно сразу сказать, что условие Куранта тут нарушено. Число Куранта > 1. В результате возник отрицательный коэффициент (-0.12) там, где (при выполнении условия Куранта) должны были быть только положительные.
Это приводит к нарушению устойчивости схемы и делает её бесполезной чуть менее, чем полностью. По результату это тоже видно, как только значения начинают превышать максимум краевых значений 1.0: при устойчивой схеме такое немыслимо Думаю, дальше там и отрицательные значения повылезают. Ошибки округления растут экспоненциально, и довольно быстро забивают всё решение. Нет смысла говорить о порядке точности в таких условиях.

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение28.05.2018, 01:25 
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 по таким же рассуждениям тоже появится отрицательный член.

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение28.05.2018, 13:27 
Аватара пользователя
Да, я наврал в предыдущем сообщении, поспешил :oops: Хотел спинным мозгом обойтись, а тут надобно головным думать :-)

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

Соответственно, и норму невязки придётся считать не как максимум $([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$-я степень шага при уменьшении шагов с сохранением условия Куранта.

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение29.05.2018, 16:57 
Благодарю, вчера еще увидел, но уснул на 15 часов :oops:

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

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

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение29.05.2018, 18:11 
Аватара пользователя
a3ro в сообщении #1315917 писал(а):
Не будет ли в данном случае халтурой, если я буду просто каждый элемент численных решений сравнивать с 0 и 1 и брать наименьшую разность, чтобы подставить в фон Неймана?
Не знаю. Боюсь, не получится.
a3ro в сообщении #1315917 писал(а):
Моё непонимание в том, что тут 2 функции, а дифференцируются по 1 переменной, а не наоборот.
Когда одна переменная, а функций две, это гораздо проще. Пара обычных дифуров, не в частных производных. Можно обойтись без разностных схем, а методами вроде Рунге-Кутта. Только проблема, что неизвестных больше, чем уравнений :-) Ну то есть одного уравнения Эйлера мало. Или вам давление дано, и вы по нему можете вычислить скорость, или наоборот. Или ещё что-то.
Когда переменных две или больше, тут уже будут уравнения в частных производных, и схема понадобится. Только та же проблема остаётся: не хватает одного этого уравнения, нужно ещё что-то, либо давление, либо скорость, либо ещё какое-то уравнение. Нужно поработать над постановкой, а потом уже задумываться над методами решения.
И потом, тут вы перескакиваете от модельной математической задачи к более сложной задаче гидродинамики, не слишком ли быстро?

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение29.05.2018, 19:18 
Цитата:
Не знаю. Боюсь, не получится.

Да просто не хочу "вручную" забивать точные данные для каждой точки в сетках, где число куранта не 0, чтобы сравнивать (хочу обойтись спинным мозгом :D)
Цитата:
Только проблема, что неизвестных больше, чем уравнений :-) Ну то есть одного уравнения Эйлера мало. Или вам давление дано, и вы по нему можете вычислить скорость, или наоборот. Или ещё что-то.Нужно поработать над постановкой, а потом уже задумываться над методами решения.

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

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

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

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

 
 
 
 Re: Начально-краевая задача для уравнения переноса.
Сообщение29.05.2018, 22:15 
Аватара пользователя
По разностной схеме для задачи гидродинамики я бы порекомендовал создать новую тему.
От заголовка этой темы она отошла.

 
 
 [ Сообщений: 16 ]  На страницу 1, 2  След.


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group