Нужна помощь при расчете по разностной схеме

[Zhenia]

Столкнулся вот с какой проблемой

Начну с постановки. Мне необходимо решить следующую задачу одномерной гидродинамики (Лагранжевы массовые координаты). Система уравнений (адиабатическое приближение)

$$ \frac{du}{dt} + \frac{dg}{ds} = 0, \quad \frac{d}{dt}\frac{1}{\rho} = \frac{du}{ds}, \quad \frac{dx}{dt} = u, \quad p = c\rho^2, \quad g = p + \omega, \quad \omega = - \nu \rho \frac{du}{ds} $$

задана в области $ 0 < s < M, \quad 0 < t < T$ при дополнительных условиях

$$ u(0,t) = 0, \quad u(M,t) = 0, \quad u(s,0) = u_0(s), \quad \rho(s,0) = 1. $$

Здесь $u(s,t)$ - скорость (обычная, не массовая), $x(s,t)$ - координата, $p(s,t)$ - давление, $\rho(s,t)$ - плотность, $\nu$коэффициент вязкости, $c$ - скорость звука в среде, $t$ - вермя. Лагранжева массовая координата задается соотношением

$$ s = \int_0^x \rho(\xi,t)\,d\xi $$

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

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

Итак, дле решения данной задачи я использую конечно-разностный подход. Из книжки Самарского и Попова [Разностные методы решения задач газовой динамики, М.:Наука 1992, стр. 200] списывается полностью консервативная неявная разностная схема. Реализуя данную схему, из системы дифференциальных уравнений получаем системы нелинейных алгебраических уравнений, которые решаем итерационным методом Ньютона. В методе Ньютона на каждой итерации решаем систему уже линейных алгебраических уравнений с трехдиагональной матрицей. Решаем, естественно, методом прогонки.

Возьмем для примера нулевые начальные условия для скорости. При счете ошибка (точное решение здесь будет u(s,t) = 0) некоторое время растет линейно (как и обещает схема - первый порядок аппроксимации по времени). Но, начиная с некоторого момента, скорость роста ошибки вдруг становится экспоненциальной и вся схема разваливается (просто, за несколько временных шагов происходит переполнение разрядной сетки). При всем при этом в том же Самарском и Попове есть доказательство абсолютной устойчивости данной схемы.

Пытаясь отыскать причину такого странного поведения, я проверил следующие места в своих расчетах

1. Корректность работы прогонки. Условия устойчивости и существования решения, получаемого прогонкой, при расчетах выполняются.
2. Закон сохранения вещества. Обеспечивается при формулировке абсолютно консервативной схемы. Численно выполняется с высокой точностью.
3. Условие сходимости итераций Ньютона. Тут самое интересное. На каждом временном шаге я проверяю условие сходимости (в равномерной норме):

$$ q = ||1/\rho^{k+1}|| \, ||\rho^{k+1} - \rho^{k}|| $$

здесь $\rho$ - сеточная функция плотности, k - номер итерации. Для сходимости величина q должна быть меньше единицы. В моих расчетах на первом шаге по времени значение q очень хорошее, но со временем оно 'портится' - растет. Экспоненциальный рост ошибки начинается именно после того, как q перевалит за единицу.
4. Что меня удручает больше всего. Уменьшая шаг по времени хотелось бы надеяться на то, что схема будет считать точнее. Однако на деле чем меньше я делаю шаг по времени, тем быстрее начинается резкий рост ошибки.

Бьюсь над задачей уже несколько месяцев, но никак не могу найти корень проблемы. Подозреваю, что либо я где-то я наврал в расчетной программе, либо здесь проблема с накоплением ошибки, которую я упускаю из вида. Может быть, необходимо округлять промежуточные результаты или же вести счет с переменным шагом по времени?

Дорогие коллеги, может, кто-нибудь сталкивался в своей практике с подобными проблемами и поделится опытом, где нужно искать ошибку? Возможно, есть литература в которой описаны различные подводные камни, ожидающие вычислителей? Буду признателен за любую помощь!

Всем спасибо!

[photon]

Zhenia писал(а):

...или же вести счет с переменным шагом по времени?

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

[G^a]

Добречко!

С ходу могу рекомендовать:

1) Почитать про некорректные, жёсткие задачи. И поняв, откуда берётся проблема ухудшения результата с уменьшением шага по времени - почитать об операциях регуляризации.
2) Поднять порядок схемы.

2) Попробовать применить метод конечных элементов - он сложнее, но зачастую более устойчив и корректен.

[photon]

G^a писал(а):

2) Попробовать применить метод конечных элементов - он сложнее, но зачастую более устойчив и корректен.

или воспользоваться готовыми пакетами, например FemLAB

[Zhenia]

Дорогие друзья!

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

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

G^a писал(а):

1) Почитать про некорректные, жёсткие задачи. И поняв, откуда берётся проблема ухудшения результата с уменьшением шага по времени - почитать об операциях регуляризации.
2) Поднять порядок схемы.

2) Попробовать применить метод конечных элементов - он сложнее, но зачастую более устойчив и корректен.

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

Всех благодарю!

[Zhenia]

Господа!

Не подскажете ли источники литературы, посвященные решению жестких задач, возникающих при численном моделировании (разностные схемы для УрЧП).

Всех благодарю!