Об аппроксимации начальных данных

В пояснениях к схеме «крест» приведена простейшая аппроксимация условия приводящая к погрешности , в то время как во внутренних узлах сетки погрешность аппроксимации для некоторых из рассматриваемых здесь схем является величиной второго порядка как по t, так и по h.

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

Представим значение решения в точках первого слоя по времени в виде ряда Тейлора по степеням t:

Замечаем, что из дифференциального уравнения следует

Таким образом,

Отсюда получаем расчетную формулу для данных на первом слое по времени:

Последнее соотношение, переписанное в виде

очевидно, аппроксимирует условие со вторым порядком точности.

О последовательности вычислений

Сначала из соотношений для начальных данных схемы «крест» находятся значения: (m = 0, …, M) и (m = 1, 2, …, M – 1). Недостающие компоненты сеточной функции на первом слое и определяются из условий на границах. Затем для p = 1, 2, …, P – 1 из соотношения (10.2) отыскиваются (m = 0, 1, …, M). Для неявных схем необходимо решать систему линейных уравнений с трехдиагональной матрицей.

Сведение задачи (10.1) к задаче для системы двух уравнений первого порядка

Замечание. В этом пункте демонстрируется, как с помощью искусственного приема рассматриваемую задачу можно свести к другой — известной как задача для уравнений акустики. Последние — ничто иное, как уравнения переноса. Таким образом, этот раздел представляет собой связующее звено между данной лабораторной работой и работой «Численное решение дифференциальных уравнений в частных производных гиперболического типа. Уравнение переноса».

Задача (10.1) эквивалентна задаче:

где функция, стоящая в правой части первого уравнения

(10.6)

Обоснование

Введем в рассмотрение новую функцию связав ее с уравнением и полагая Тогда волновое уравнение можно записать в виде

Интегрируем последнее соотношение по t, получаем:

После выполнения интегрирования

т. е. где

Таким образом, задача (10.1) сведена к задаче (10.6) для двух дифференциальных уравнений первого порядка.

Дополнительные замечания

Легко проверить непосредственно, что дифференциальные уравнения (10.6) можно записать в виде:

где — инварианты Римана.

Особенность последних уравнений состоит в том, что в каждом из них дифференцируется лишь одна неизвестная функция (q или r). Отметим, что возможность записи исходной системы типа (10.6) в виде системы для инвариантов (функций, для которых уравнения разделяются так же, как и здесь) является признаком гиперболичности системы.

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

(10.7)

Уравнения для инвариантов, к которым мы перешли, являются известными уравнениями переноса (см. лабораторную работу «Численное решение дифференциальных уравнений в частных производных гиперболического типа. Уравнение переноса»). Последнее обстоятельство позволит распространить разностные схемы для уравнения переноса на наш случай.

Таким образом, вместо исходной задачи (10.1) можно решать задачу (10.6) или задачу (10.7) для инвариантов. В последнем случае исходная искомая функция u вычисляется через q и r в нужных точках по формуле

Варианты разностных схем для задачи (10.6)

Предварительные замечания. Как известно, неявные схемы практически всегда устойчивы. Однако, применительно к задаче (10.6) использование их связано с определенными трудностями, так как расчет данных на очередном слое по времени требует решения линейной системы уравнений (более общего вида, нежели система с трехдиагональной матрицей). Поэтому в рамках данной работы мы сосредоточим внимание на способах построения более простых (с точки зрения реализации) явных разностных схем для задачи (10.6). С ними, однако, применительно к задаче (10.6) снова не все ясно. (Попробуйте по произвольному «явному» шаблону построить устойчивую схему!). Поэтому сначала разберемся как строятся подходящие явные схемы для задачи (10.7), сформулированной для инвариантов q и r.

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

СХИ1 — схема для инвариантов (первого порядка точности).Эта схема является естественным обобщением на случай системы для инвариантов схем типа «явный угол» для уравнения переноса. С учетом характеристических направлений (направлений переноса инвариантов) первое из уравнений (10.9) аппроксимируется по шаблону

шаблон для второго уравнения

Разностные уравнения для внутренних узлов сетки:

p = 1, 2, …, P – 1; m = 0, 1, …, M – 1;

(10.8)

p = 1, 2, …, P – 1; m = 1, 2, …, M.

Аппроксимация краевых условий:

p = 1, 2, …, P – 1. (10.9)

Начальные данные:

m = 0, 1, …, M. (10.10)

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

Последовательность вычислений для схемы СХИ1.Из соотношений (10.10) находятся q и r во всех точках начального слоя. Затем в рамках стандартного программного цикла для p = 0, 1, 2, …, P – 1 осуществляется расчет данных на (p + 1)-м слое по времени. При этом:

1) из соотношений (10.8) находятся для m = 0, …, M – 1 и для m = 1, 2, …, M;

2) из (10.9) с использованием найденных значений и отыскиваются и

СХИ2 — схема для инвариантов (второго порядка точности).В этой схеме используется шаблон

Замечание. Данная схема является обобщением на случай системы (10.7) явной четырехточечной схемы второго порядка для уравнения переноса.

Разностные уравнения для внутренних узлов сетки:

p = 1, 2, …, P – 1; m = 1, 2, …, M – 1. (10.11)

Аппроксимация краевых условий:

p = 1, 2, …, P – 1. (10.12)

Начальные данные:

m = 0, 1, …, M. (10.13)

В отличие от схемы СХИ1 данная система соотношений пока не замкнута. Подходящий способ замыкания данной схемы рассматривается ниже. При выбранном там способе замыкания схема имеет второй порядок точности и устойчива при выполнении условия Куранта:

Дополнительные соотношения схемы СХИ2. Заметим, что дифференциальное уравнение для инварианта q — ничто иное, как характеристическое соотношение

выполняющееся вдоль характеристик которые представляют собой семейство прямых

В частности, это характеристическое соотношение связывает значение q в точках левой границы со значениями q внутри расчетной области (на — характеристике, выпущенной из рассматриваемой точки левой границы). Поэтому логично для получения одного недостающего уравнения аппроксимировать с желаемым порядком точности дифференциальное уравнение для q из (10.7) в ячейках сетки, примыкающих к левой границе. Аналогичным образом второе дополнительное соотношение получается аппроксимацией уравнения для r из (10.7) по ячейке, примыкающей к правой границе.

Замечание. В схеме СХИ1 разностные уравнения для q при m = 0 и для r при m = M, можно рассматривать как аппроксимации первого порядка характеристических соотношений: для q вблизи левой границы и для r вблизи правой.

В этом случае разностные уравнения (10.11) аппроксимируют соответствующие дифференциальные уравнения для инвариантов со вторым порядком точности. Поэтому желательно, чтобы необходимые замыкающие соотношения обеспечивали ту же точность.

Приведем соответствующий пример явной аппроксимации дифференциального уравнения для q вблизи левой границы.

На рис. 10 изображен фрагмент сеточной области вблизи левой границы Точка А является точкой пересечения — характеристики, выпущенной из узла B с координатами с предыдущим слоем по времени.

Интегрируя характеристическое соотношение

по отрезку AB, получим

Заменим интеграл, например, по формуле трапеций. Имеем:

Значение определим по известным значениям q в узлах p-го слоя по времени с помощью интерполяционной формулы:

здесь d = at — расстояние точки A от левой границы.

Подставляя выражение для в предыдущую формулу, получаем одно из недостающих дополнительных соотношений:

Другое недостающее уравнение получим, аппроксимируя аналогичным образом характеристическое соотношение для инварианта r вблизи правой границы:

Расчет по схеме СХИ2. Из соотношений (10.13) находятся q и r во всех точках начального слоя. Затем в рамках стандартного программного цикла (для p = 0, 1, 2, …, P – 1) осуществляется расчет данных на -м слое по времени. При этом:

1) из уравнений (10.11) находится а также для m = 1, 2, …, M – 1;

2) из уравнений (10.12) с использованием найденных из дополнительных соотношений и отыскиваются и

Теперь рассмотрим примеры явных разностных схем непосредственно для задачи (10.6).

Схема СХ1. В этой схеме используется шаблон

Разностные уравнения для внутренних узлов сетки:

p = 1, 2, …, P – 1; m = 1, 2, …, M – 1. (10.14)

Аппроксимация краевых условий:

p = 1, 2, …, P – 1. (10.15)

Начальные данные:

m = 0, 1, …, M. (10.16)

Как и в случае схемы СХИ2 данная система соотношений пока не замкнута. Подходящий способ замыкания данной схемы обсуждается ниже. При разумном способе замыкания схема имеет первый порядок точности и устойчива, если выполнено условие Куранта:

Дополнительные соотношения для схемы СХ1.Недостающие для вычисления значений в точках левой и правой границы соотношения можно получить, аппроксимируя, например, дифференциальное уравнение для v по тому или иному шаблону вблизи левой и правой границ. Возможны следующие варианты.

A) При m = 0: шаблон «явный правый уголок», уравнение

При m = M: шаблон «явный левый уголок», уравнение

Дополнительные соотношения:

и

B) При m = 0: шаблон «неявный левый уголок», уравнение

При m = M: шаблон «явный правый уголок», уравнение

Дополнительные соотношения:

Замечание 1. Возникают естественные вопросы: почему выбрано дифференциальное уравнение для при конструировании дополнительных соотношений? Почему выбраны те или иные шаблоны? Точного ответа на первый вопрос нет. Что касается шаблона, то тут руководствуемся тем, чтобы не ухудшить аппроксимацию и устойчивость схемы в целом. Вариант A), по-видимому, непригоден, так как аппроксимация уравнений (10.7) по таким шаблонам, вообще говоря, неустойчива.

Замечание 2. Данная схема является «близкой родственницей» схемы СХИ1 для инвариантов. В самом деле, уравнения группы (10.14) этой схемы являются линейной комбинацией (суммой и разностью) соответствующих уравнений схемы СХИ1. Тем самым проясняется вопрос о происхождении данной схемы.

Учитывая отмеченное обстоятельство, можно сделать вывод, что естественным способом замыкания данной схемы является привлечение в качестве дополнительных соотношений уравнений для q (при m = 0) и r (при m = M) из группы уравнений (10.14) схемы СХИ1. В этом случае, очевидно, рассматриваемая схема алгебраически эквивалентна схеме СХИ1 для инвариантов.

Схема СХ2.В этой схеме используется шаблон

Разностные уравнения для внутренних узлов сетки представляют собой линейные комбинации (полусумму и полуразность) уравнений для q и rсхемы СХИ2 для инвариантов:

p = 1, 2, …, P – 1; m = 1, 2, …, M – 1. (10.17)

Краевые условия (10.15) и начальные данные (10.16) записываются так же как в схеме СХ1.

Как и в случае схемы СХИ2 для инвариантов, данная система соотношений пока не замкнута. Подходящий способ замыкания обсуждается ниже. При соответствующем способе замыкания схема имеет второй порядок точности и устойчива, если выполнено условие Куранта:

Дополнительные соотношения для схемы СХ2.Недостающие для вычисления значений v в точках левой и правой границы соотношения можно получить, как и для схемы СХ1, аппроксимируя, например, дифференциальное уравнение для v по тому или иному шаблону вблизи левой и правой границ.

Способы, использованные для замыкания схемы СХ1, приводят к погрешности аппроксимации первого порядка, в то время, как уравнения для внутренних узлов данной схемы обеспечивают второй порядок аппроксимации.

В качестве подходящего способа получения дополнительных уравнений, не портящих второй порядок точности схемы в целом, применим аппроксимацию уравнения для v системы (10.6) вблизи левой и правой границы по прямоугольному шаблону.

При m = 0 имеем

Отсюда получаем дополнительную расчетную формулу, из которой находится

Аналогичным образом, при m = M получаем еще одну недостающую формулу (для ).

Применительно к данной схеме справедливы замечания, сделанные при обсуждении способов замыкания схемы СХ1. Естественным способом представляется привлечение характеристических соотношений. В частности, если использовать в качестве дополнительных уравнений аппроксимацию характеристических соотношений вблизи границ со вторым порядком точности, то данная схема будет алгебраически эквивалентна СХ1.

Расчет по схемам СХ1 и СХ2. Из соотношений (10.16) находятся u и во всех точках начального слоя. Затем осуществляется расчет данных на (p + 1)-м слое по времени для значений p = 0, …, P – 1. При этом:

1) из уравнений (10.14) находятся величины и для m = 1, 2, …, M – 1;

2) из (10.15) с использованием дополнительных соотношений отыскиваются и