При классификации уравнений с частными производными (2.1) отмечалось, что уравнения первого порядка называются также уравнениями переноса. Это объясняется тем, что такие уравнения описывают процессы переноса частиц в средах, распространения возмущений и т.п.
В общем случае уравнения переноса могут иметь значительно более сложный вид (например, интегродифференциальное уравнение Больцмана в кинетической теории газов). Однако здесь мы ограничимся линейным уравнением с частными производными первого порядка. Его решение представляет интерес не только с практической точки зрения; в еще большей степени это уравнение полезно при разработке и исследовании разностных схем.
Будем считать, что искомая функция Uзависит от времени tи одной пространственной переменной х. Тогда линейное уравнение переноса может быть записано в виде
(2.23)
Здесь а — скорость переноса, которую будем считать постоянной и положительной. Это соответствует переносу (распространению возмущений) слева направо в положительном направлении оси х. Правая часть F(x, t) характеризует наличие поглощения (или, наоборот, источников) энергии, частиц и т.п. в зависимости от того, какой физический процесс описывается уравнением переноса.
Характеристики уравнения (2.23) определяются соотношениями х — at = С = const. При постоянном а они являются прямыми линиями, которые в данном случае (а > 0) наклонены вправо (рис. 2.5).
Рис. 2.5. Область решения
Расчетная область при решении уравнения (2.23) может быть как бесконечной, так и ограниченной. В первом случае, задавая начальное условие при t = 0:
(2.24)
получаем задачу Коши для полуплоскости На практикеобычно приходится решать уравнение переноса в некоторой ограниченной области (например, в прямоугольнике ; см. рис. 2.5). Начальное условие (2.24) в этом случае задается на отрезке l1; граничное условие нужно задать при х = 0, т.е. на отрезке l2, поскольку при а > 0 возмущения распространяются вправо. Это условие запишем в виде
(2.25)
Таким образом, задача состоит в решении уравнения (2.23) с начальным и граничным условиями (2.24) и (2.25) в ограниченной области G:
Убедиться в том, что данная задача поставлена правильно (корректно) можно, проанализировав решение уравнения (2.23), которое при F(x, t) = 0 имеет вид
(2.26)
где Н — произвольная дифференцируемая функция. В этом легко убедиться, подставляя (2.26) в уравнение (2.23). Решение (2.26) называется бегущей волной (со скоростью а). Это решение постоянно вдоль каждой характеристики: при х — at = С искомая функция U = Н(х — at) = Н(С) постоянна. Таким образом, начальные и граничные условия переносятся вдоль характеристик, поэтому они должны задаваться на отрезках l1и l2 расчетной области G(см. рис. 2.5).
Можно также построить аналитическое решение задачи Коши для неоднородного уравнения (2.23). Заметим лишь, что решение этой задачи меняется вдоль характеристики, а не является постоянным.
Рассмотрим разностные схемы для решения задачи (2.23) — (2.25). Построим в области Gравномерную прямоугольную сетку с помощью прямых xi = ih (i =0,1. I) и . Вместо функций U(x,t), F(x,t), Ф(х) и будем рассматривать сеточные функции, значения которых в узлах (xi, tj) соответственно равны и . Для построения разностной схемы необходимо выбрать шаблон. Примем его в виде правого нижнего уголка(рис. 2.6). При этом входящие в уравнение (2.23) производные аппроксимируются конечно-разностными соотношениями с использованием односторонних разностей:
(2.27)
Рис. 2.6. Правый нижний уголок
Решая это разностное уравнение относительно единственного неизвестного значения на (j + 1)-ом слое, получаем следующую разностную схему:
(2.28)
Полученная схема явная, поскольку значения сеточной функции в каждом узле верхнего слоя выражаются явно с помощью соотношений (2.28) через ранее найденные ее значения на предыдущем слое.
Для начала счета по схеме (2.28), т.е. для вычисления сеточной функции на первом слое, необходимы ее значения на слое j= 0. Они определяются начальным условием (2.24), которое записываем для сеточной функции:
(2.29)
Граничное условие (2.25) также записывается в сеточном виде:
(2.30)
Таким образом, решение исходной дифференциальной задачи (2.23) — (2.25) сводится к решению разностной задачи (2.28) – (2.30). Найденные значения сеточной функции принимаются в качестве значений искомой функции и в узлах сетки.
Алгоритм решения исходной задачи (2.23) — (2.25) с применением рассмотренной разностной схемы достаточно прост. На рис. 2.7 представлена его структурограмма. В соответствии с этим алгоритмом в памяти компьютера хранится весь двумерный массив , и он целиком выводится на печать по окончании счета. С целью экономии памяти (и если эти результаты не понадобятся для дальнейшей обработки) можно воспользоваться тем, что схема двухслойная, и хранить лишь значения сеточной функции на двух соседних слоях . Рекомендуем читателю соответственным образом модифицировать представленный алгоритм и построить новую структурограмму.
Рис. 2.7. Алгоритм решения линейного уравнения переноса
Укажем теперь некоторые свойства данной разностной схемы. Она аппроксимирует исходную задачу с первым порядком, т.е. невязка имеет порядок O(h+τ). Схема условно устойчива; условие устойчивости имеет вид
(2.31)
Эти свойства схемы установлены в предположении, что решение U(x, t), начальное и граничное значения Ф(х) и дважды непрерывно дифференцируемы, а правая часть F(x, t) имеет непрерывные первые производные.
Поскольку схема (2.28) устойчива и аппроксимирует исходную задачу, то в соответствии с приведенной в разд. 2.1 теоремой сеточное решение сходится к точному с первым порядком при . Отметим, что при а 0 эта схема не сходится.
Граничное условие для уравнения переноса (2.23) при а 0). Такая аппроксимация называется противопотоковой и широко используется при численном решении уравнений переноса.
При построении явной разностной схемы (2.28) производная ¶U/¶х аппроксимировалась с помощью значений сеточной функции на j-ом слое; в результате получилось разностное уравнение (2.27), в котором использовано значение сеточной функции лишь в одном узле верхнего слоя. Если производную¶U/¶х аппроксимировать на (j + 1)-ом слое (шаблон изображен на рис. 2.9), то получится неявная схема. Разностное уравнение примет вид
(2.34)
Рис. 2.9. Правый верхний уголок
Разрешая это уравнение относительно , приходим к следующей разностной схеме:
(2.35)
Это двухслойная трехточечная схема первого порядка точности. Она безусловно устойчива (при а > 0). Хотя формально данная разностная схема строилась как неявная, практическая организация счета по ней проводится так же, как и для явных схем.
Действительно, в правую часть уравнения (2.35) входит значение на (j+1)-ом слое, которое при вычислении уже найдено. При расчете значение берется из граничного условия (2.30). По объему вычислений и логике программы (см. рис. 2.7) схема (2.35) аналогична схеме (2.28), однако безусловная устойчивость делает ее более удобной, поскольку исключается ограничение на величину шага.
Схему (2.28) можно применять для решения задачи Коши в неограниченной области, поскольку граничное условие (2.30) в этой схеме можно не использовать.
Рис. 2.10. Прямоугольник
Рассмотрим еще одну разностную схему, которую построим на симметричном прямоугольном шаблоне (рис. 2.10). Производная по tздесь аппроксимируется в виде полусуммы отношений односторонних конечных разностей в (i — 1)-м и i-м узлах, а производная по x — в виде полусуммы конечно-разностных соотношений на j—ми (j + 1)-ом слоях. Правую часть вычисляют в центре ячейки, хотя возможны и другие способы ее вычисления (например, в виде некоторой комбинации ее значений в узлах). В результате указанных аппроксимаций получим разностное уравнение в виде
(2.36)
Данная двухслойная четырехточечная схема также формально построена как неявная. Однако из (2.36) можно выразить неизвестное значение через остальные, которые предполагаются известными:
(2.37)
Построенная схема имеет второй порядок точности. Она устойчива на достаточно гладких решениях.
Схема (2.37) получена для случая а > 0. Аналогичную ей схему при а 0, а2 > 0 — скорости переноса вдоль осей х, у, (2.39) — начальное условие при t= 0; (2.40) — граничные условия при х =0, y= 0.
В трехмерной области (х, у, t) построим разностную сетку, ячейки которой имеют форму прямоугольного параллелепипеда. Для этого проведем координатные плоскости через точки деления осей х, у, t:
Значение сеточной функции в узле (i, j, k), с помощью которой аппроксимируются значения , обозначим через . Построим безусловно устойчивую разностную схему первого порядка точности, аналогичную схеме (2.35). Шаблон изображен на рис. 2.11, где выделена одна ячейка разностной сетки. Сплошными линиями соединены узлы шаблона. Нижний слой (нижнее основание параллелепипеда) имеет номер k, верхний k+ 1.
Рис. 2.11. Шаблон для двумерного уравнения
По аналогии с (2.34) запишем разностное уравнение, аппроксимирующее дифференциальное уравнение (2.38):
Разрешим это уравнение относительно значения сеточной функции в узле :
(2.41)
Вычислительный алгоритм этой схемы аналогичен алгоритму одномерной схемы (2.35). Здесь также счет производится по слоям k= 1,2. К. При k= 0 используется начальное условие (2.39), которое нужно переписать в разностном виде:
(2.42)
На каждом слое последовательно вычисляют значения сеточной функции в узлах. При этом последовательность перехода от узла к узлу может быть различной: двигаются параллельно либо оси х, либо оси у. Во втором случае последовательность вычисляемых значений следующая:
На рис. 2.12 показана нумерация узлов, соответствующая данной последовательности вычислений на каждом временном слое. Точками отмечены расчетные узлы сетки, крестиками — граничные узлы, в которых значения сеточной функции задаются граничными условиями (2.40). Эти условия обходимо записать в сеточном виде:
. (2.43)
Рис. 2.12. Последовательность вычислений
При этом значения в угловой точке (х = 0, у = 0) в данной разностной схеме не используются.
Алгоритм решения смешанной задачи (2.38 – 2.40) для двумерного уравнения переноса по схеме (2.41) с учетом сеточных начального и граничных условий (2.42) и (2.43) представлен на рис. 2.13. При этом некоторые блоки (вычисление начальных значений uij, значений на границе пересылка ) даны схематически, хотя каждый из них представляет циклический алгоритм.
Рис. 2.13. Алгоритм решения двумерного уравнения переноса
В данном алгоритме предусмотрено хранение в памяти машины не полного трехмерного массива искомых значений , а лишь значений на двух слоях: — нижний слой, — верхний слой (искомые значения). Введен счетчик выдачи l, решение выдается через каждые Lслоев; при L = 1 происходит выдача результатов на каждом слое. Блок «Вычисление » вычисляет искомое значение по формуле, которая в принятых в структурограмме обозначениях имеет вид
Видео:27. Уравнения переносаСкачать
Дипломная работа: Разностные схемы для уравнения переноса на неравномерных сетках
Название: Разностные схемы для уравнения переноса на неравномерных сетках Раздел: Рефераты по информатике, программированию Тип: дипломная работа Добавлен 06:44:11 13 ноября 2009 Похожие работы Просмотров: 1664 Комментариев: 21 Оценило: 3 человек Средний балл: 5 Оценка: неизвестно Скачать | |||||||||||
Рис.2,а. Трехточечный шаблон | Рис.2,б. Трехточечный шаблон | Рис.2,в. Трехточечный шаблон | Рис.2,г. Четырехточечный шаблон |
Составим разностные схемы ко всем четырем шаблонам на рис.2.
(7а)
, (7б)
, (7в)
. (7г)
Во всех четырех схемах правая часть выбиралась в центре ячейки. Возможен и другой способ аппроксимации правой части.
Все четыре разностные схемы (7а) — (7г), по существу, являются явными. Во всех схемах значение явно выражается через . Решение на нулевом слое известно из начального условия, т.е. . Для вычисления решения на следующем слое из граничного условия находим , это позволяет найти , далее вычисляется и т.д. Таким образом находится решение на первом слое, аналогично находится решение на втором слое и т.д. Именно в связи с тем, что решение вычисляется слой за слоем слева направо, схемы (7а) — (7г) называются схемами бегущего счета.
Алгоритмы бегущего счета обеспечивают существование и единственность решений при любых . Поэтому для доказательства сходимости остается разобраться с аппроксимацией и устойчивостью разностных схем. Поскольку граничное условие воспроизводится точно, постольку исследование устойчивости по нему не требуется.
Разностная схема (7а). Исследуем погрешность аппроксимации схемы (7а). Для этого разложим решение и правую часть в окрестности точки (tm,xn) в ряд Тейлора, считая непрерывность всех требуемых производных:
,
,
.
Учитывая эти разложения, находим невязку схемы (7а):
т.е. схема (7а) имеет аппроксимацию первого порядка в норме .
Устойчивость исследуем с помощью принципа максимума, формулировка и доказательство которого приведены в лекции №9. Критерий равномерной устойчивости по начальным данным (формула (64) в лекции №9 при C = 0) дает следующее ограничение:
,
которое сводится к так называемому условию Куранта
Согласно (8), разностная схема (7а) является условно устойчивой© в норме .
Методом разделения переменных можно доказать необходимость условия (8) для обеспечения устойчивости. Подставим в схему (7а) следующие величины:
,
тогда множитель роста гармоники
.
Условие устойчивости обеспечивается, когда
. (9)
Выполнение неравенства (9) при произвольном q обеспечено, когда r £ 1, т.е. при выполнении условия Куранта. При нарушении условия Куранта, т.е. при r > 1 неравенство (9) не выполняется при всех q, а только при некоторых. Так, при r >> 1 неравенство (9) перепишется в виде: cos qh £ ½, т.е. амплитуды некоторых гармоник растут при переходе со слоя на слой и схема неустойчива по начальным данным.
Устойчивость по правой части согласно формуле (65) лекции №9 обеспечивается при k = 1 в норме , когда верно условие Куранта.
В итоге схема (7а) при выполнении условия Куранта сходится в с первым порядком точности.
В качестве примера рассмотрим численное решение задачи
(10)
Задача (10) имеет следующее аналитическое решение:
(11)
На листинге_№1 приведен код программы численного решения задачи (10) по разностной схеме (7а). На рис.3,а приведено трехмерное изображение решения u(t,x) при выполнении условия Куранта, а на рис.3,б приведено решение при нарушении условия Куранта. Видно, появление неустойчивости в решении при нарушении условия (8).
%Программа численного решения уравнения
%очищаем рабочее пространство
%определяем параметр скорости переноса c,
%а также отрезок времени интегрирования T и
%диапазон изменения пространственной
%определяем шаг по пространству
%рассматривается два варианта расчета
%при tau=h/c (условие Куранта выполняется) и
%при tau=1.12*h/c (условие Куранта нарушено)
%определяем сетки по времени и по пространству
%определяем начальное значение u(0,x)=x^3/(12c^2)
%организуем расчет по разностной схеме (7а)
%определяем левое граничное значение
%рисуем численное решение уравнения переноса u(t,x)
surf(ti,xi,y); [xi ti]=meshgrid(x,t);
%рисуем численное решение уравнения переноса u(t,x)
Рис.3,а. Численное решение уравнения (10) по разностной схеме (7а) при выполнении условия Куранта | Рис.3,б. Численное решение уравнения (10) по разностной схеме (7а) с нарушением условия Куранта (8) |
Сравним теперь численное решение задачи (10) и аналитическое решение (11). На листинге_№2 приведен код соответствующей программы. В программе считается, что t = 0.5h/c и варьируется шаг по пространству. На рис.4 приведен итог работы кода программы листинга_№2 в виде кривой зависимости отношения ошибки численного решения к шагу сетки const(h) = в зависимости от шага сетки h. Из условия аппроксимации разностной схемой (7а) исходного уравнения (3) с порядком O(t + h) следует, что величина const(h) должна стремиться к некоторой константе по мере того, как h ® 0. Такая тенденция видна на рис.4.
%Программа численного решения уравнения
%сравнение его с аналитическим решением
%определяем параметр скорости переноса c,
%а также отрезок времени интегрирования T и
%диапазон изменения пространственной
%определяем количество делений шага h пополам
%делим шаг h пополам
%определяем шаг по времени, который считается
%пропорциональным шагу по пространству
%определяем сетки по времени и по пространству
%определяем начальное значение u(0,x)=x^3/(12c^2)
%организуем расчет по разностной схеме (7а)
%определяем левое граничное значение
%определяем ошибку численного решения в норме C
%и делим ее на шаг сетки h
%рисуем зависимость предстепенной константы от
%функция, возвращающая аналитическое решение
Разностная схема (7б) исследуется аналогично. Для исследования аппроксимации разложение в ряд Тейлора удобно проводить в окрестности узла (xn —1,tm + t). Для дважды непрерывно дифференцируемого решения данная схема при выполнении условия устойчивости
обеспечивает сходимость со скоростью O(t + h).
Разностная схема (7в) безусловно устойчива и на дважды непрерывно дифференцируемых решениях сходится к точному решению со скоростью O(t + h).
Разностная схема (7г) симметричная и следует ожидать, что порядок ее аппроксимации выше, чем в предыдущих членах. Для оценки порядка аппроксимации разложение в ряд Тейлора удобно провести в окрестности центра ячейки . После проведения соответствующих выкладок, можно найти оценку невязки:
. (13)
Тем самым схема (7г) имеет второй порядок аппроксимации, когда решения имеют непрерывные производные вплоть до третьей.
Рис.4. Зависимость предстепенной константы в оценке ошибки
численного решения от шага сетки
Устойчивость разностной схемы (7г) исследуем с помощью метода разделения переменных. Подставляя в (7г)
,
найдем значение коэффициента роста Фурье-гармоники при переходе со слоя на слой:
. (14)
Из оценки (14) следует, что для любой гармоники и при любых соотношениях шагов. Это означает, что схема (7г) безусловно и равномерно устойчива по начальным данным в норме .
Исследуем разностную схему (7г) на предмет сходимости в двух нормах: и . На листинге_№3 приведен код программы для изучения сходимости схемы (7г) на примере численного решения задачи (10) и сравнения полученного решения с аналитическим решением (11). В программе вычисляются зависимости предстепенных констант const(h) для двух норм от шага сетки h, при этом считается, что t = 0.5h/c. Согласно теоретическим оценками, предстепенная константа const(h) = должна выходить на некоторое постоянное значение при h ® 0.
%Программа численного решения уравнения
%сравнение его с аналитическим решением
%определяем параметр скорости переноса c,
%а также отрезок времени интегрирования T и
%диапазон изменения пространственной
%определяем количество делений шага h пополам
%делим шаг h пополам
%определяем шаг по времени, который считается
%пропорциональным шагу по пространству
%определяем сетки по времени и по пространству
%определяем начальное значение u(0,x)=x^3/(12c^2)
%организуем расчет по разностной схеме (7г)
🔍 Видео
Лекция 291. Схема ускоренного переносаСкачать
Лекция №1.1 Явная и неявная схемы для уравнения теплопроводностиСкачать
Лукьяненко Д. В. - Численные методы - Лекция 18Скачать
Лукьяненко Д. В. - Численные методы - Лекция 19Скачать
6-2. Метод сетокСкачать
№6. Уравнения в частных производных. Уравнения переноса, мелкой воды.Скачать
Кобельков Г. М. - Численные методы. Часть 2. Семинары - Лекция 6Скачать
Вычислительная математика 17 Теория разностных схемСкачать
Кобельков Г. М. - Численные методы. Часть 2 - Решение уравнений в частных производныхСкачать
Лукьяненко Д. В. - Численные методы - Лекция 22Скачать
Решение уравнений. Как переносить слагаемые из одной части уравнения в другую. Математика 6 классСкачать
Построение логических схемСкачать
Теоретические основы электротехники 24. Расчёт схем с помощью теории графов, топологических матриц.Скачать
Как решать уравнения по схеме ГорнераСкачать
Решение задачи теплопроводности (Явная разностная схема)Скачать
Сеточные методы решения дифференциальных уравнений в частных производных.Скачать
Лекция 8. Булева интерпретация релейных схемСкачать
Вычислительная математика 18 Численные методы решения уравнений в частных производныхСкачать