Часовой пояс: UTC + 3 часа [ Летнее время ]

Численные методы решения уравнений в частных производных
Онлайн-сервисы Нахождение НОД и НОК Разложение числа на простые множители Сравнения по модулю Операции над множествами Операции над векторами Разложение вектора по базису. Доказательство, что векторы образуют базис Чертёж треугольника по координатам вершин Решение треугольника Решение Пирамиды Построение Пирамиды по координатам вершин Чертёж многоугольника по координатам вершин Решение систем методом Крамера и Матричным Онлайн построение графика кривой 2-го порядка Определение вида кривой или поверхности 2-го порядка по инвариантам МНК и регрессионный анализ Онлайн + графики Онлайн число, сумма и дата прописью Алгоритмы JavaScript Алгоритмы поиска Алгоритмы сортировки Уникальные элементы массива Объединение, пересечение и разность массивов НОД и НОК Операции над матрицами Дата прописью Введение в анализ Функции: понятие, определение, графики Непрерывность функции Исследование функции и построение графика Теория множеств Множества: понятие, определение, примеры Точечные множества Замкнутые и открытые множества Мера множества Группы, кольца, поля в математике Поле комплексных чисел Кольцо многочленов Основная теорема алгебры и ее следствия Математическая логика Алгебра высказываний Аксиоматика и логические рассуждения Методы доказательств теорем Алгебра высказываний и операции над ними Формулы алгебры высказываний Тавтологии алгебры высказываний Логическая равносильность формул Нормальные формы для формул высказываний Логическое следование формул Приложение алгебры высказываний для теорем Дедуктивные и индуктивные умозаключения Решение логических задач Принцип полной дизъюнкции Булевы функции Множества, отношения и функции в логике Булевы функции от одного и двух аргументов Булевы функции от n аргументов Системы булевых функций Применение булевых функций к релейно-контактным схемам Релейно-контактные схемы в ЭВМ Практическое применение булевых функций Теория формального Формализованное исчисление высказываний Полнота и другие свойства формализованного исчисления высказываний Независимость системы аксиом формализованного исчисления высказываний Логика предикатов Логика предикатов Логические операции над предикатами Кванторные операции над предикатами Формулы логики предикатов Тавтологии логики предикатов Преобразования формул и следование их предикатов Проблемы разрешения для общезначимости и выполнимости формул Применение логики предикатов в математике Строение математических теорем Аристотелева силлогистика и методы рассуждений Принцип полной дизъюнкции в предикатной форме Метод полной математической индукции Необходимые и достаточные условия Логика предикатов и алгебра множеств Формализованное исчисление предикатов Неформальные и формаль- ные аксиоматические теории Неформальные аксиоматические теории Свойства аксиоматических теорий Формальные аксиоматические теории Формализация теории аристотелевых силлогизмов Свойства формализованного исчисления предикатов Формальные теории первого порядка Формализация математической теории Теория алгоритмов Интуитивное представление об алгоритмах Машины Тьюринга и тезис Рекурсивные функции Нормальные алгоритмы Маркова Разрешимость и перечислимость множеств Неразрешимые алгоритмические проблемы Теорема Гёделя о неполноте формальной арифметики Математическая логика и компьютеры Математическая логика и языки программирования Применение компьютеров для доказательства теорем математической логики Математическая логика и логическое программирование Математическая логика и информатика Математическая логика и искусственный интеллект Дискретная математика Множества и отношения Теория множеств: понятия и определения Операции над множествами Кортеж и декартово произведение множеств Соответствия и бинарные отношения на множествах Операции над соответствиями на множествах Семейства множеств Специальные свойства бинарных отношений Отношения эквивалентности на множестве Упорядоченные множества Теорема о неподвижной точке Мощность множества Парадокс Рассела Метод характеристических функций Группы и кольца Алгебраические структуры и операции Группоиды, полугруппы, группы Кольца, тела, поля Области целостности в теории колец Модули и линейные пространства Подгруппы и подкольца Теорема Лагранжа о порядке конечной группы Гомоморфизмы групп и нормальные делители Гомоморфизмы и изоморфизмы колец Алгебра кватернионов Полукольца и булевы алгебры Полукольца: определение, аксиомы, примеры Замкнутые полукольца Полукольца и системы линейных уравнений Булевы алгебры и полукольца Решетки и полурешетки Алгебраические системы Алгебраические системы: модели и алгебры Подсистемы алгебраических систем Конгруэнции и фактор-системы Гомоморфизмы алгебраических систем Прямые произведения алгебраических систем Конечные булевы алгебры Многосортные алгебры Теория графов Теория графов: основные понятия и определения Способы представления графов Неориентированные и ориентированные деревья Остовное дерево и алгоритм Краскала Методы систематического обхода вершин графа Алгоритмы поиска в глубину и ширину в графах Задача о путях во взвешенных ориентированных графах Изоморфизм, гомоморфизм и автоморфизм графов Топологическая сортировка вершин графа Элементы цикломатики в теории графов Булева алгебра и функции Булевы функции и булев куб Таблицы булевых функций и булев оператор Равенство булевых функций. Фиктивные переменные Формулы и суперпозиции булевых функций Дизъюнктивные и конъюнктивные нормальные формы Построение минимальных ДНФ Теорема Поста и классы Критерий Поста Схемы из функциональных элементов Конечные автоматы и регулярные языки Конечные автоматы и регулярные языки Алфавит, слово, язык в программировании Порождающие грамматики (грамматики Хомского) Классификация грамматик и языков Регулярные языки и регулярные выражения Конечные автоматы Допустимость языка конечным автоматом Теорема Клини Детерминизация конечных автоматов Минимизация конечных автоматов Лемма о разрастании для регулярных языков Обоснование алгоритма детерминизации автоматов Конечные автоматы с выходом Морфизмы и конечные подстановки Машины Тьюринга Контекстно-свободные языки Контекстно-свободные языки и грамматики Приведенная форма КС-грамматики Лемма о разрастании для КС-языков Магазинные автоматы (автомат с магазинной памятью) Алгоритм построения МП-автомата по КС-грамматике Алгоритм построения КС-грамматики по МП-автомату Алгебраические свойства КС-языков Основное свойство суперпозиции КС-языков Пересечение контекстно-свободных языков Методы синтаксического анализа КС-языков Восходящий синтаксический анализ и LR(k)-грамматики Семантика формальных языков Принцип индукции по неподвижной точке Графовое представление МП-автоматов Интегральное исчисление Неопределённый и определённый Неопределенный и определенный интегралы Свойства интегралов Интегрирование по частям Интегрирование методом замены переменной Интегрирование различных рациональных функций Интегрирование различных иррациональных функций Интегрирование различных тригонометрических функций Определенный интеграл и его основные свойства Необходимое и достаточное условие интегрируемости Теоремы существования первообразной Свойства определенных интегралов Несобственные интегралы Интегральное определение логарифмической функции Приложения интегралов Вычисление площадей плоских фигур Площади фигур в различных координатах Вычисление объемов тел с помощью интегралов Объём тела вращения Вычисление длин дуг кривых Формулы длины дуги регулярной кривой Кривизна плоской кривой Площадь поверхности вращения тела Интегралы в физике Статические моменты и координаты центра тяжести Теоремы Гульдина–Паппа Вычисление моментов инерции Другие приложения интегралов в физике Основные интегралы Интеграл Ньютона-Лейбница Интеграл Римана Интеграл Лебега Вариационное исчисление Примеры вариационных задач Дифференциальное уравнение Эйлера Функционалы, зависящие от нескольких функций Задача о минимуме кратного интеграла Финансовый анализ Анализ эффективности Критерии и показатели эффективности предприятия Методы анализа эффективности деятельности Факторный анализ прибыли от операционной деятельности Анализ безубыточности предприятия Операционный рычаг и эффект финансового рычага Анализ и оценка состава, структуры и динамики доходов и расходов Анализ рентабельности и резервов устойчивого роста капитала Анализ распределения прибыли предприятия Анализ и оценка чувствительности показателей эффективности Анализ устойчивости Финансовая устойчивость и долгосрочная платежеспособность Характеристика типов финансовой устойчивости Рыночная активность Финансовый анализ рыночной активности Методика анализа рыночной активности Анализ и оценка дивидендного дохода на одну акцию Инвестиционная деятельность Инвестиции: экономическая сущность и классификация Государственное регулирование инвестиционной деятельности Источники финансовых ресурсов на капитальные вложения Инвестиции в основные фонды Оценка состояния основных фондов Амортизация основных фондов Капитальное строительство в инвестиционном процессе Планирование инвестиций в форме капитальных вложений Экономическая эффективность инвестиций Финансирование капитальных вложений Кредитование капитальных вложений Кредитоспособность Финансирование и кредитование затрат Финансирование и кредитование инвестиционной деятельности потребительской кооперации Финансирование и кредитование капитальных вложений потребительской кооперации Инвестиционное строительное проектирование Анализ инвестиций Инвестиции и инвестиционная деятельность предприятия Задачи финансового анализа инвестиций предприятия Учет фактора времени в инвестиционной деятельности Аннуитет и финансовая рента в инвестициях Учет фактора инфляции при инвестировании Оценка фактора риска инвестиционного проекта Методы оценки эффективности инвестиций Показатели эффективности инвестиционного проекта Стоимость компании Концепция построения международных стандартов финансовой отчетности (МСФО) Экономическое содержание международных стандартов финансовой отчётности Цели и принципы оценки стоимости акций и активов компании Оценка акций и активов предприятия по справедливой стоимости Методы оценки справедливой стоимости акций предприятия Затратный подход к оценки стоимости компаний и акций Сравнительный подход к оценки стоимости предприятий и акций Доходный подход к оценке стоимости компании и акций Выбор ставки дисконтирования при инвестировании в акции Метод капитализации прибыли Сравнение подходов к оценке стоимости компаний и пакетов акций Форвардные контракты Форвардный контракт и цена Форвардная цена акции на бирже Цена форвардного контракта инвестора Форвардная цена акции с учетом величины дивиденда Форвардная цена акции с учетом ставки дивиденда Форвардная цена валюты на рынке форекс Форвардный валютный курс и инфляция на рынке Форвардная цена товара и спотовый рынок Форвардная цена при различии ставок по кредитам и депозитам Синтетический форвардный контракт на акции и валюту Теория вероятностей Основные понятия теории вероятностей Зависимые и независимые случайные события Повторные независимые испытания Формула Бернулли Одномерные случайные величины Многомерные случайные величины Функции случайных величин Законы распределения целочисленных случайных величин Законы распределения непрерывных случайных величин Предельные теоремы теории вероятностей Закон больших чисел и предельные теоремы Вероятностные закономерности Математическая статистика Элементы математической статистики Выборочный метод Оценки параметров генеральной совокупности Статистические гипотезы Критерии согласия Теоретические и эмпирические частоты Теория очередей (СМО) Определение системы массового обслуживания Уравнения Колмогорова Предельные вероятности состояний Определение СМО с отказами Определение СМО с ожиданием (очередью) Аналитическая геометрия Векторная алгебра Метрические понятия и аксиомы геометрии Равенство и подобие геометрических фигур Бинарные отношения Вектор, его направление и длина Линейные операции над векторами Линейная зависимость и независимость векторов Отношение коллинеарных векторов Проекции векторов на прямую и на плоскость Угол между векторами Ортогональные проекции векторов Координата вектора на прямой и базис Координаты вектора на плоскости и базис Координаты вектора в пространстве и базис Операции над векторами в координатной форме Ортогональный и ортонормированный базисы Cкалярное произведение векторов и его свойства Выражение скалярного произведения через координаты векторов Векторное произведение векторов и его свойства Смешанное произведение векторов и его свойства Ориентированные площади и объемы Двойное векторное произведение и его свойства Применение векторов в задачах на аффинные свойства фигур Применение произведений векторов при решении геометрических задач Применение векторной алгебры в механике Системы координат Прямоугольные координаты Преобразования прямоугольных координат Полярная система координат Цилиндрическая система координат Сферические координаты Аффинные координаты Аффинные преобразования координат Аффинные преобразования плоскости Примеры аффинных преобразований плоскости Аффинные преобразования пространства Многомерное координатное пространство Линейные и аффинные подпространства Скалярное произведение n-мерных векторов Преобразования систем координат Геометрия на плоскости Алгебраические линии на плоскости Общие уравнения геометрических мест точек Алгебраические уравнения линий на плоскости Уравнения прямой, проходящей через точку перпендикулярно вектору Уравнения прямой, проходящей через точку коллинеарно вектору Уравнения прямой, проходящей через две точки Уравнения прямой с угловым коэффициентом Взаимное расположение прямых Примеры задач с прямыми на плоскости Системы неравенств с двумя неизвестными Системы линейных уравнений с двумя неизвестными Линии 2-го порядка Канонические уравнения линий второго порядка Порядок приведения уравнения линии к каноническому виду Эллипс Гипербола Парабола Квадратичные неравенства с двумя неизвестными Применение линий 1-го и 2-го порядков в задачах на экстремум функций Инварианты линий Классификация линий 2-го порядка по инвариантам Приведение уравнения линии к каноническому виду по инвариантам Геометрия в пространстве Способы задания ГМТ в пространстве Алгебраические уравнения поверхностей Уравнения плоскости, проходящей через точку перпендикулярно вектору Уравнения плоскости, компланарной двум неколлинеарным векторам Уравнения плоскости, проходящей через три точки Взаимное расположение плоскостей Типовые задачи с плоскостями Уравнения прямых в пространстве Взаимное расположение прямых в пространстве Типовые задачи с прямыми в пространстве Поверхности 2-го порядка Канонические уравнения поверхностей Порядок приведения уравнения поверхности к каноническому виду Поверхности второго порядка Эллипсоиды Гиперболоиды Конусы Параболоиды Применение поверхностей 1-го и 2-го порядков в задачах на экстремум функций Инварианты поверхностей Классификация поверхностей 2-го порядка по инвариантам Квадратичные неравенства с тремя неизвестными Приведение уравнения поверхности к канониче-скому виду по инвариантам Линейная алгебра Матрицы и операции Линейные операции над матрицами Умножение матриц Возведение матриц в степень Многочлены от матриц Транспонирование и сопряжение матриц Блочные матрицы Произведение и сумма матриц Кронекера Метод Гаусса приведения матрицы к ступенчатому виду Элементарные преобразования матриц Определители Определители матриц и их основные свойства Формула полного разложения определителя Формула Лапласа полного разложения определителя Определитель произведения матриц Методы вычисления определителей Ранг матрицы Линейная зависимость и линейная независимость строк (столбцов) матрицы Ранг матрицы и базисный минор матрицы Методы вычисления ранга матрицы Ранг системы столбцов (строк) Обратная матрица Обратные матрицы и их свойства Ортогональные и унитарные матрицы Способы нахождения обратной матрицы Матричные уравнения Односторонние обратные матрицы Скелетное разложение матрицы Полуобратная матрица Псевдообратная матрица Системы уравнений Системы линейных алгебраических уравнений Метод Гаусса решения систем линейных уравнений Структура общего решения системы уравнений Решение систем с помощью полуобратных матриц Псевдорешения системы линейных уравнений Функциональные матрицы Функциональные матрицы скалярного аргумента Производные матриц по векторному аргументу Линейные и квадратичные формы и их преобразования Приведение форм к каноническому виду Закон инерции вещественных квадратичных форм Знакоопределенность форм вещественных квадратичных Формы и исследование функций на экстремум Многочленные матрицы Многочленные матрицы (лямбда-матрицы) Операции над лямбда-матрицами Простые преобразования многочленных матриц Инвариантные множители многочленной матрицы Функции от матриц Собственные векторы и значения матрицы Подобие числовых матриц Характеристический многочлен матрицы Минимальный многочлен матрицы Теорема Гамильтона-Кэли Жорданова форма матрицы Приведение матрицы к жордановой форме Многочлены от матриц Применение многочленов от матриц Функции от матриц Линейные пространства Линейные пространства: определение и примеры Линейная зависимость и независимость n-мерных векторов Размерность и базис линейного пространства Преобразования координат в линейном пространстве Изоморфизм линейных пространств Подпространства Подпространства линейного пространства Пересечение и сумма подпространств Способы описания подпространств Нахождение дополнения и суммы подпространств Нахождение пересечения подпространств Линейные отображения Линейные многообразия Линейные отображения Матрица линейного отображения Ядро и образ линейного отображения Линейные операторы Линейные операторы (преобразования) Инвариантные подпространства Собственные векторы и значения оператора Свойства собственных векторов операторов Канонический вид линейного оператора Методика приведения линейного преобразования к каноническому виду Евклидовы пространства Евклидовы пространства Ортогональные векторы евклидова пространства Ортогональный базис евклидова пространства Ортонормированный базис евклидова пространства Ортогональные дополнения в евклидовом пространстве Задача о перпендикуляре Матрица и определитель Грама и его свойства Линейные преобразования евклидовых пространств Канонический вид ортогонального оператора евклидова пространства Сопряженные операторы евклидова пространства Самосопряженные операторы евклидова пространства Приведение квадратичной формы к главным осям Унитарные пространства и их линейные преобразования Комплексный анализ Комплексные числа Комплексные числа в алгебраической форме Комплексные числа в тригонометрической и показательной формах Множества на комплексной плоскости Последовательности и ряды комплексных чисел Комплексные функции Функции комплексного переменного. Предел, непрерывность и производная Элементарные функции комплексного переменного Дифференцирование функций комплексного переменного Аналитические функции и их свойства Конформные отображения и их свойства Интегрирование функций комплексного переменного Функциональные ряды в комплексной области Функциональные ряды и последовательности Степенные ряды и их свойства Разложение функций в степенные ряды Нули аналитических функций Ряд Лорана и разложение функций по целым степеням Особые точки, Вычеты Изолированные особые точки функций и полюсы Вычеты и их применение Вычисление интегралов с помощью вычетов Вычеты и расположение нулей многочлена Операционное исчисление Преобразование Лапласа и его свойства Решение ДУ операционным методом Анализ выходных процессов линейных стационарных систем Z-преобразование и его свойства Дифференциальные уравнения ДУ первого порядка Основные понятия и определения ДУ Метод изоклин для ДУ 1-го порядка Метод последовательных приближений ДУ с разделяющимися переменными Однородные ДУ Линейные ДУ 1-го порядка Дифференциальное уравнение Бернулли ДУ в полных дифференциалах Интегрирующий множитель ДУ, не разрешенные относительно производной Дифференциальное уравнение Риккати Составление ДУ семейств линий Задачи на траектории Особые решения ДУ ДУ высших порядков Понятия и определения ДУ высших порядков ДУ, допускающие понижение порядка Линейная независимость функций Определители Вронского и Грама Однородные и неоднородные дифференциальные уравнения Задача Коши и Уравнение Эйлера Линейные ДУ с переменными коэффициентами Метод Лагранжа решения ДУ Краевые задачи для ДУ высших порядков Разложение решения ДУ в степенной ряд Разложение решения ДУ в обобщенный степенной ряд Нахождение периодических решений ДУ Асимптотическое интегрирование ДУ Системы ДУ Системы ДУ: понятия и определения Сведение системы ДУ к одному уравнению Нахождение интегрируемых комбинаций Интегрирование однородных линейных систем ДУ Методы интегрирования неоднородных систем ДУ Преобразование Лапласа и решение ДУ и систем Теория устойчивости Устойчивость решений ДУ по Ляпунову Простейшие типы точек покоя Метод функций Ляпунова Устойчивость решений ДУ по первому приближению Критерии устойчивости Рауса–Гурвица и Михайлова ДУ с малым параметром при производной Численные методы Методы алгебры Численные методы линейной алгебры Численные методы решения СЛАУ Итерационный метод Шульца обратной матрицы Методы решения задач о собственных значениях и векторах матрицы Методы решения нелинейных уравнений Методы решения систем нелинейных уравнений Методы теории приближений Методы приближения сеточных функций Методы функциональной интерполяции Методы интегрально-дифференциальной интерполяции Методы интегрального сглаживания Методы интерполяции и сглаживания сплайнами Методы численного дифференцирования и интегрирования Методы численного дифференцирования Методы численного интегрирования Методы решения обыкновенных ДУ Численные методы решения задачи Коши Разностные схемы для решения задачи Коши Составные схемы для решения задачи Коши Экстраполяционные методы решения задачи Коши Непрерывно-дискретные методы решения задачи Коши Численные методы решения краевых задач Методы решения ДУ в частных производных Численные методы решения уравнений математической физики с двумя переменными Принципы построения разностных схем для уравнений в частных производных Разностные схемы решения уравнений в частных производных 1-го порядка Разностные схемы решения уравнений в частных производных 2-го порядка Численные методы решения уравнений в частных производных Численные методы решения уравнений математической физики с тремя переменными	Численные методы решения уравнений в частных производных Метод прямых Основная идея метода состоит в сведении проблемы решения уравнений в частных производных к решению системы обыкновенных дифференциальных уравнений. Он используется в основном для решения уравнений эллиптического и параболического типов. В методе прямых, в отличие от метода сеток, не все дифференциальные операторы, входящие в исходную дифференциальную задачу, заменяются разностными. Обычно для аппроксимации выбирается оператор, соответствующий пространственной переменной в задаче теплопроводности или соответствующий неравноправному направлению в задачах с эллиптическими уравнениями. Тем самым дифференциальная задача заменяется дифференциально-разностной с меньшим числом непрерывных переменных (понижается размерность дифференциальной задачи). В связи с этим данный метод можно рассматривать как предельный случай метода сеток, когда шаги сетки по определенным независимым переменным стремятся к нулю. Основную идею конструирования метода прямых рассмотрим на примере нахождениия решения $u(x,y)$ уравнения: $\begin{gathered}A(x,y)\frac{\partial^2 u(x,y)}{\partial x^2}+ 2B(x,y)\frac{\partial^2 u(x,y)}{\partial y\,\partial x}+ C(x,y) \frac{\partial^2 u(x,y)}{\partial y^2}+ D(x,y)\frac{\partial u(x,y)}{\partial x}\,+\\[2pt] +\,E(x,y)\frac{\partial u(x,y)}{\partial y}+ F(x,y)u(x,y)= G(x,y),\quad (x,y)\in \Omega= (0,L)\times (0,M), \end{gathered}$ (8.59) с краевыми условиями $\begin{aligned}& (x,0)= \varphi(x), & & u(x,M)= \psi(x),\\ & (0,y)= \Phi(y), & & u(L,y)= \Psi(y) \end{aligned}$ (8.60) на границах области прямоугольной формы (рис. 8.15), где $\varphi(x),\psi(x), \Phi(x), \Psi(x)$ — заданные функции. Пусть выполнено условие $B^2-AC<0$ эллиптичности уравнения для всех $(x,y)\in \Omega$ . Отсюда следует, что $A(x,y)\ne 0,~ C(x,y)\ne 0$ . Для решения задачи выберем на отрезке $[0,M]$ точки $y_{j}= j\,h,~ 0 \leqslant j \leqslant J,~ J\,h=M$ , и проведем через них прямые, параллельные оси $Ox$ . На каждой прямой дифференциальное уравнение в частных производных (8.59) приближенно заменим обыкновенным дифференциальным уравнением для функций $u(x,y_{j})$ . Для этого аппроксимируем частные производные по $y\colon$ $\begin{gathered}\left.{\frac{\partial \widehat{u}(x,y)}{\partial y}}\right\|_{y=y_{j}}= \frac{\widehat{u}(x,y_{j+1})-\widehat{u}(x,y_{j-1})}{2h}\,,\\[2pt] \left.{\frac{\partial^2 \widehat{u}(x,y)}{\partial y^2}}\right\|_{y=y_{j}}= \frac{\widehat{u}(x,y_{j+1})-2\,\widehat{u}(x,y_{j})+ \widehat{u}(x,y_{j-1})}{h^2}\,,\\[2pt] \left.{\frac{\partial^2 \widehat{u}(x,y)}{\partial y\, \partial x}}\right\|_{y=y_{j}}= \frac{\widehat{u}_{x}(x,y_{j+1})-\widehat{u}_{x}(x,y_{j-1})}{2h}\,. \end{gathered}$ (8.61) Введем обозначения $\begin{gathered}\widehat{u}(x,y_{j})= \widehat{u}_{j}(x),\qquad \frac{\partial \widehat{u}(x,y_{j})}{\partial x}= \widehat{u}\,'_{j}(x),\qquad \frac{\partial^2 \widehat{u}(x,y_{j}) }{\partial x^2}= \widehat{u}\,''_{j}(x),\\[4pt] A(x,y_{j})= A_{j}(x),\qquad B(x,y_{j})= B_{j}(x),\qquad C(x,y_{j})= C_{j}(x)\quad \text{etc.}\end{gathered}$ (8.62) На прямых, описываемых уравнениями $y_{j}=\text{const},~ j=1,\ldots,J-1$ , получаем $\begin{gathered} A_{j}(x)\widehat{u}\,''_{j}(x)+ B_{j}(x)\frac{\widehat{u}\,'_{j+1}(x)-\widehat{u}\,'_{j-1}(x)}{h}+ C_{j}(x) \frac{\widehat{u}_{j+1}(x)-2\,\widehat{u}_{j}(x)+ \widehat{u}_{j-1}(x)}{h^2}+ D_{j}(x) \widehat{u}\,'_{j}(x)+\\[2pt] E_{j}(x) \frac{\widehat{u}_{j+1}(x)-\widehat{u}_{j-1}(x)}{2h}+ F_{j}(x) \widehat{u}_{j}(x)= G_{j}(x),\quad j=1,\ldots,J-1;~ 0<x<L. \end{gathered}$ Здесь $\widehat{u}_{1}(x),\ldots, \widehat{u}_{J-1}(x)$ — приближенное решение задачи на прямых $y=y_1=\text{const},\ldots, y=y_{J-1}= \text{const},$ $A_{j}(x),\ldots, G_{j}(x)$ — коэффициенты уравнения (8.59) на этих прямых. В силу краевых условий находим $\begin{gathered}\widehat{u}_{0}(x)= u(x,0)= \varphi(x),\qquad \widehat{u}_{J}(x)= u(x,M)= \psi(x),\\[4pt] \widehat{u}_{j}(0)= \Phi(y_{j})= \Phi(j\,h),\quad \widehat{u}_{j}(L)= \Psi(y_{j})= \Psi(j\,h),\quad j=1,\ldots,J-1. \end{gathered}$ (8.63) Таким образом, от линейного дифференциального уравнения в частных производных (8.59) и краевыми условиями (8.60) переходим к дифференциально-разностной краевой задаче (8.62),(8.63) для системы из $(n-1)$ обыкновенных дифференциальных уравнений (8.62) относительно $(n-1)$ неизвестных функций $\widehat{u}_{1}(x),\ldots, \widehat{u}_{J-1}(x)$ . В результате ее решения можно найти законы изменения искомой функции $u(x,y)$ вдоль прямых $y_{j}=\text{const},$ $j=0,1,\ldots,J$ , то есть $\widehat{u}(x,y_{j}),~ j=0,1,\ldots,J$ . Значения функции в промежуточных точках по у может быть найдено с помощью интерполяции. Система (8.62),(8.63) называется системой уравнений метода прямых. Замечания 1. В зависимости от выбора аппроксимационных формул для производных по у можно получить различные системы метода прямых, отличающиеся порядком аппроксимации. 2. Метод прямых удобнее применять в том случае, когда коэффициенты в уравнении (8.59) не зависят от $x$ . Тогда система (8.62) будет системой с постоянными коэффициентами. 3. Если расчетная область отлична от прямоугольника, при реализации метода возникают трудности, связанные с тем, что функции $u_{j}(x)$ определены на различных множествах. 4. При решении нестационарных уравнений дискретизация проводится по пространственным переменным, а время $t$ остается непрерывной переменной. 5. На практике система дифференциальных уравнений (8.62) из-за большого числа $J$ должна решаться численно, поэтому потребуется вводить шаг не только по $y$ . 6. Если изменение искомой функции внутри множества $\Omega$ в направлении $y$ не велико, то можно ограничиться одной или двумя прямыми. Пример 8.9. Составить уравнения метода прямых для решения одномерной задачи теплопроводности (где $\psi(x)$ — заданная функция) $\begin{aligned}& u_{t}=a^2\cdot u_{xx},\quad 0<x<1,\quad 0<t<1,\\ & u(x,0)= \psi(x),\quad 0 \leqslant x \leqslant 1,\\ & u(0,t)= u(1,t)=0,\quad 0<t \leqslant 1. \end{aligned}$ Решение, формула (8.64) Для решения задачи выберем на отрезке $[0,1]$ точки $x_{i}= i\,h,~ i=0,1,\ldots, I;~ I\,h=1$ и проведем через них прямые, параллельные оси $Ot$ (рис. 8.16). На каждой прямой приближенно заменим уравнение теплопроводности обыкновенным дифференциальным уравнением для функции $u(x_{i},t)$ . Для этого заменим частную производную $u_{xx}$ по формуле $\left.{\frac{\partial^2 \widehat{u}}{\partial x^2}}\right\|_{x=x_{i}}= \frac{\widehat{u}(x_{i-1},t)-2\,\widehat{u}(x_{i},t)+ \widehat{u}(x_{i+1},t)}{h^2}$ и введем обозначения $\widehat{u}(x_{i},t)= \widehat{u}_{i}(t),~ \frac{\partial \widehat{u}(x_{i}, t)}{\partial t}= \widehat{u}\,'_{i}(t)$ . Тогда на прямых $x_{i}=\text{const},~ i=1,\ldots,I-1$ , получаем систему и $\begin{gathered}\widehat{u}\,'_{i}(t)= a^2\cdot \frac{\widehat{u}(x_{i-1},t)-2\,\widehat{u}(x_{i},t)+ \widehat{u}(x_{i+1},t)}{h^2}\,,\quad 1 \leqslant i \leqslant I-1,\\[4pt] \begin{aligned}& \widehat{u}_{0}(t)= u(0,t)=0, & & 0<t \leqslant 1,\\ & \widehat{u}_{I}(t)= u(1,t)=0, & & 0<t \leqslant 1 \end{aligned} \end{gathered}$ (8.64) с начальными условиями $\widehat{u}_{i}(0)= \psi(x_{i})= \psi(i\,h),~ i=1,2,\ldots,I-1$ . Полученную систему дифференциальных уравнений можно записать в векторно-матричной форме: $\frac{d}{dt}\! \begin{pmatrix}\widehat{u}_1\\ \vdots\\ \vdots\\ \widehat{u}_{I-1} \end{pmatrix}= \frac{a^2}{h^2}\! \begin{pmatrix}-1&1&\cdots& 0&0\\ 1&-2&1&\cdots&0\\ \vdots& \vdots& \vdots& \vdots& \vdots\\ \vdots& \vdots& 1&-2&1\\ 0&0&\cdots& 1&-2 \end{pmatrix}\!\cdot\! \begin{pmatrix}\widehat{u}_1\\ \vdots\\ \vdots\\ \widehat{u}_{I-1}\end{pmatrix}\!.$ Задача имеет аналитическое решение $\widehat{u}(t)= H^{-1}\cdot\! \begin{pmatrix}\exp \dfrac{\lambda_1a^2t}{h^2}&0&0\\ 0&\ddots&0\\ 0&0& \exp \dfrac{\lambda_{I-1}a^2t}{h^2} \end{pmatrix}\!\cdot H\,u(0),$ где $u(0)= \begin{pmatrix}\psi(h)\\ \psi(2h)\\ \vdots\\ \psi \bigl((I-1)h\bigr)\end{pmatrix}\!,~ \lambda_{i}=-\left(2\cos \frac{\pi\,i}{n}+ 2\right)\!,~ i=1,\ldots,I-1$ — собственные значения матрицы системы, элементы матрицы $H$ (собственные векторы матрицы $A$ ) равны $H_{ik}= \sin \frac{\pi\,ik}{n}\,,$ $i=1,\ldots,I-1,$ $k=1,\ldots,I-1$ . Заметим, что при численном решении системы обыкновенных дифференциальных уравнений при увеличении числа $I$ (улучшении точности) коэффициент жесткости $S$ системы увеличивается. При больших значениях $n$ справедливо $S= \frac{\max_{i}\|\lambda_{i}\|}{\min_{i}\|\lambda_{i}\|}\cong \frac{4I^2}{\pi^2}$ . Метод характеристик Данный метод применяется при численном решении уравнений математической физики, описывающих различные процессы в газовой динамике, теории упругости, прочности и других прикладных областях. При его использовании в процессе решения задачи имеется возможность выделять в расчетных областях особые поверхности или линии, отделяющие области с одними свойствами от областей с другими свойствами. Рассмотрим основные идеи метода на примере задачи Коши для квазилинейной системы дифференциальных уравнений m-го порядка с двумя независимыми переменными $x,\,t\colon$ $\begin{aligned}&\begin{aligned}&a_{11} \frac{\partial u_{1}(x,t)}{\partial t}+ \ldots+a_{1m} \frac{\partial u_{m}(x,t)}{\partial t}+ b_{11} \frac{\partial u_{1}(x,t)}{\partial x}+\ldots+b_{1m} \frac{\partial u_{m}(x,t)}{\partial x}= f_1(x,t,u),\\[-2pt] &\quad\vdots\\[-4pt] &a_{m1} \frac{\partial u_{1}(x,t)}{\partial t}+ \ldots+a_{mm} \frac{\partial u_{m}(x,t)}{\partial t}+ b_{m1} \frac{\partial u_{1}(x,t)}{\partial x}+\ldots+b_{mm} \frac{\partial u_{m}(x,t)}{\partial x}= f_m(x,t,u),\end{aligned}\\[4pt] &\begin{aligned}& u_1(x,0)= \varphi_1(x), & & x\in \mathbb{R},\\[-2pt] & \quad\vdots\\[-4pt] & u_m(x,0)= \varphi_m(x), & & x\in \mathbb{R}, \end{aligned} \end{aligned}$ где $a_{ij}=a_{ij}(x,t,u),~ b_{ij}=b_{ij}(x,t,u),~ i=1,\ldots,m;~ j=1,\ldots,m;~ D=R\times (0,T),~ T>0$ — заданный момент времени. Перепишем систему в векторно-матричном виде $\begin{gathered}A(x,t,u)\frac{\partial u(x,t)}{\partial t}+ B(x,t,u)\frac{\partial u_{1}(x,t)}{\partial x}= f(x,t,u),\quad (x,t)\in D,\\[2pt] u(x,0)= \varphi(x),\quad x\in \mathbb{R}, \end{gathered}$ (8.65) где $A(x,t,u)= \bigl(a_{ij}(x,t,u)\bigr),~ B(x,t,u)= \bigl(b_{ij}(x,t,u)\bigr)$ — матрицы размеров $(m\times m),~ f(x,t,u)= \begin{pmatrix}f_1(x,t,u)\\\vdots\\ f_m(x,t,u) \end{pmatrix}$ — вектор-функция, $u(x,t)= \begin{pmatrix}u_1(x,t)\\ \vdots\\ u_m(x,t) \end{pmatrix}$ — искомая вектор-функция правых частей, $\varphi(x)= \begin{pmatrix} \varphi_1(x)\\ \vdots\\ \varphi_m(x) \end{pmatrix}$ — заданная вектор-функция. Пусть известно, что система (8.65) имеет гладкое решение в области $D$ и на некоторой кривой $\gamma$ оно задано и равно $\Bigl.{u(x,t)}\Bigr\|_{\gamma}$ (рис. 8.17). Рассмотрим задачу о нахождении решения в окрестности кривой $\gamma$ . Выберем на кривой $\gamma$ произвольную точку $(x_0,t_0)$ , а вектор бесконечно малого смещения вдоль этой кривой из точки $(x_0,t_0)$ обозначим $(dx,dt)$ . Присоединим к уравнениям (8.65) выражение для полного дифференциала (для упрощения записи опустим аргументы): $A\,\frac{\partial u}{\partial t}+ B\,\frac{\partial u}{\partial x}=f,\qquad du= \frac{\partial u}{\partial x}\,dx+ \frac{\partial u}{\partial t}\,dt\,.$ Умножим первое уравнение на $dt$ , а из второго выразим $\frac{\partial u}{\partial t}\,dt\colon$ $A\,\frac{\partial u}{\partial t}\,dt+ B\,\frac{\partial u}{\partial x}\,dt= f\,dt,\qquad \frac{\partial u}{\partial t}\,dt= du-\frac{\partial u}{\partial x}\,dx\,.$ Подставляя второе уравнение в первое, имеем $A\,du-A\,\frac{\partial u}{\partial x}\,dx+ B\,\frac{\partial u}{\partial x}\,dt= f\,dt$ или $(B\,dt-A\,dx)\cdot \frac{\partial u}{\partial x}= f\,dt-A\,du\,.$ (8.66) Если определитель матрицы $(B\,dt-A\,dx)$ отличен от нуля в каждой точке кривой $\gamma$ , то частная производная $\frac{\partial u}{\partial x}$ и, следовательно, частная производная $\frac{\partial u}{\partial t}$ определяются однозначно. В этом случае можно найти решение в окрестности кривой $\gamma\colon\, \Bigl.{u(x,t)= u(x,t)}\Bigr\|_{\gamma}+ du$ , что соответствует разложению искомого решения по формуле Тейлора до членов первого порядка. Если определитель матрицы $(B\,dt-A\,dx)$ в каждой точке кривой $\gamma$ равен нулю, то в силу предположения о существовании решения производные $\frac{\partial u}{\partial x},\, \frac{\partial u}{\partial t}$ и, следовательно, функция $u(x,t)$ , находятся неоднозначно. В этом случае кривая у называется характеристикой. Приравнивая определитель матрицы $(B\,dt-A\,dx)$ к нулю и предполагая, что $dt\ne0$ вдоль кривой $\gamma$ , получаем соотношение $\left\|B-A\cdot \frac{\partial u}{\partial t}\right\|=0\quad \text{ili}\quad \|B-A \lambda\|=0,$ (8.67) где $\lambda(x,t,u)= \frac{\partial u}{\partial t}$ . Линии, задаваемые дифференциалами смещения $dx,\,dt$ , вдоль которых справедливо (8.67), являются характеристиками. Анализ разрешимости уравнения (8.67) позволяет ввести классификацию систем вида (8.65). Если уравнение (8.67) имеет $m$ различных действительных корней, то система (8.65) называется гиперболической; если их число меньше $m$ , то система называется параболической; если уравнение не имеет действительных корней, то — эллиптической. В дальнейшем будем рассматривать только гиперболические системы, т.е. имеющие $m$ различных характеристик: $\left(\frac{dx}{dt}\right)_{i}= \lambda_{i}(x,t,u),\quad i=1,2,\ldots,m,$ (8.68) где функции $\lambda_{i}(x,t,u)$ задают характеристические направления в каждой точке кривой $\gamma$ . Вдоль характеристик решение системы должно удовлетворять характеристическим соотношениям (соотношениям совместности), которые получаются в результате приравнивания к нулю определителя m-го порядка расширенной матрицы, т.е. $\bigl(B\,dt-A\,dx,\, f\,dt-A\,du\bigr).$ (8.69) При этом следует избегать случаев получения тривиальных тождеств. Характеристические соотношения связывают дифференциалы искомых функций с дифференциалами независимых переменных вдоль характеристик. Этот способ получения характеристических соотношений основан на применении теоремы Кронекера-Капелли: решение системы (8.66) существует тогда и только тогда, когда ранг матрицы $(B\,dt-A\,dx)$ этой системы равен рангу расширенной матрицы. Предполагается, что решение существует и определитель матрицы системы (m-го порядка) вдоль характеристик равен нулю. Поэтому должен быть равен нулю определитель расширенной матрицы, составленный из произвольных $m$ столбцов. В результате описан принцип получения характеристических уравнений (уравнений характеристик и характеристических соотношений). Их вывод составляет первый этап метода характеристик. Основные идеи первого, а также второго этапа, связанного с конструированием численной схемы, рассмотрим на примере линейной задачи Коши для одномерного волнового уравнения, описывающего распространение плоских звуковых волн в покоящейся среде. Пример 8.10. Получить соотношения метода характеристик для решения линейной задачи Коши $\frac{\partial^2 \varphi(x,t)}{\partial t^2} = \frac{\partial^2 \varphi(x,t)}{\partial x^2}\,,\quad (x,t)\in \mathbb{R}\times (0;T),\quad \varphi(x,0)=x^2+x,\quad \varphi_t(x,0)=\cos x,\quad x\in \mathbb{R}.$ Решение, формулы (8.70)-(8.77) Данная задача, в частности, описывает процесс распространения звуковых волн в спокойной среде (отсутствуют источники и стоки), следовательно, решаемое уравнение можно рассматривать как уравнение акустики. Для применения метода характеристик сначала сведем задачу к виду (8.65), т.е. к системе двух уравнений первого порядка. Для этого введем новые искомые функции $\nu(x,t)= \frac{\partial \varphi(x,t)}{\partial t}\,\qquad p(x,t)= \frac{\partial \varphi(x,t)}{\partial x}\,.$ Здесь $\nu(x,t),\,p(x,t)$ — отклонения скорости и давления от их значений в невозмущенной среде, вызванные распространением звуковых волн, $\varphi(x,t)$ — потенциал функции $\nu(x,t)$ . Тогда уравнение можно представить в форме $\frac{\partial }{\partial t}\! \left(\frac{\partial \varphi(x,t)}{\partial t}\right)= \frac{\partial }{\partial x}\! \left(\frac{\partial \varphi(x,t)}{\partial x}\right)\quad \text{ili}\quad \frac{\partial \nu(x,t)}{\partial t}-\frac{\partial p(x,t)}{\partial x}=0.$ Предполагая, что задача имеет гладкое решение (т.е. рассматриваемая функция имеет столько непрерывных производных, сколько нужно для законности приводимых преобразований), воспользуемся свойством, справедливым для дважды непрерывно дифференцируемых функций: $\frac{\partial^2 \varphi(x,t)}{\partial x\,\partial t}= \frac{\partial^2 \varphi(x,t)}{\partial t\,\partial x}\quad \text{ili}\quad \frac{\partial p(x,t)}{\partial t}= \frac{\partial \nu(x,t)}{\partial x}\,.$ В результате от исходной задачи для волнового уравнения переходим к задаче для системы двух дифференциальных уравнений с начальными условиями относительно новых искомых функций: $\left\{\begin{aligned}& \frac{\partial \nu(x,t)}{\partial t}-\frac{\partial p(x,t)}{\partial x}=0, & & \nu(x,0)=\cos x,\\ & \frac{\partial p(x,t)}{\partial t}-\frac{\partial \nu(x,t)}{\partial x}=0, & & p(x,0)=2x+1. \end{aligned}\right.$ Запишем систему в виде (8.65) $A\cdot \frac{\partial u}{\partial t}+ B\cdot \frac{\partial u}{\partial x}=0,\quad u(x,0)= \begin{pmatrix}\cos x\\ 2x+1\end{pmatrix}\!,$ (8.70) где $A=\begin{pmatrix}1&0\\0&1\end{pmatrix}\!,~ B=\begin{pmatrix}0&-1\\-1&0\end{pmatrix}\!,~ u=\begin{pmatrix}\nu\\p \end{pmatrix}\!,~ m=2$ . Для решения задачи (8.70) применим метод характеристик, содержащий два этапа. Первый этап. Нахождение характеристик и вывод характеристических соотношений. Для этого решим уравнение (8.67), записанное для (8.70): $\|B-A\cdot \lambda\|= \begin{vmatrix}-\lambda&-1\\-1&-\lambda \end{vmatrix}=0.$ Получаем $\lambda^2-1=0,~ \lambda_{1,2}=\pm1$ , т.е. два направления характеристик: ${\left(\frac{dx}{dt}\right)\!}_{1}= \lambda_1=1;\quad {\left(\frac{dx}{dt} \right)\!}_{2}= \lambda_2=-1.$ (8.71) Интегрируя, найдем общие решения полученных дифференциальных уравнений. Для первого уравнения $x=t+\text{const}$ , а для второго $x=-t+\text{const}$ . Им соответствуют два семейства характеристик (в данной задаче прямых линий). Их принято обозначать буквой $C^{+}$ с верхним индексом "+" для семейства характеристик, имеющих меньший угол с положительным направлением оси Ох по сравнению с соответствующим углом для характеристик семейства $C^{-}$ . В данной задаче характеристики семейства $C^{+}$ имеют положительный наклон, а характеристики семейства $C^{-}$ — отрицательный наклон (рис. 8.18). Условно их также называют характеристиками первого и второго семейств соответственно. Для получения характеристических соотношений запишем расширенную матрицу (в решаемой задаче $f=\begin{pmatrix}0\\0 \end{pmatrix}$ ): $\bigl(B\,dt-A\,dx,\, f\,dt-A,\,du\bigr)= \left(\underbrace{\begin{pmatrix}0&-dt\\-dt&0 \end{pmatrix}}_{B\,dt}-\underbrace{\begin{pmatrix}dx&0\\ 0&dx \end{pmatrix}}_{A\,dx},-\underbrace{\begin{pmatrix}1&0\\ 0&1 \end{pmatrix}\!\! \begin{pmatrix}d\nu\\ dp\end{pmatrix}}_{A\,du} \right)= \begin{pmatrix}-dx&-dt-d\nu\\-dt&-dx&-dp \end{pmatrix}\!.$ Приравнивая определитель матрицы, составленной из последних двух столбцов, нулю, получаем $\begin{vmatrix}-dt&-d\nu\\-dx&-dp \end{vmatrix}= dp\cdot dt-d\nu\cdot dx=0.$ Поделив на $dt$ (предполагается, что $dt\ne 0$ ), имеем $dp-d\nu\, \frac{dx}{dt}=0$ . Так как на прямых семейства $C^{+}$ выполняется равенство ${\left( \frac{dx}{dt} \right)\!}_{1}= \lambda_1=1$ , то соответствующие характеристическое соотношение имеет вид $dp-d\nu=0.$ (8.72) Аналогично на прямых семейства $C^{-}$ справедливо ${\left( \frac{dx}{dt} \right)\!}_{2}= \lambda_2=-1$ и $dp+ d\nu=0.$ (8.73) Обыкновенные дифференциальные уравнения (8.71)–(8.73) образуют систему уравнений метода характеристик. Полученные характеристические соотношения (8.72), (8.73) в данной задаче могут быть проинтегрированы. В результате получаются конечные соотношения, называемые инвариантами Римана. Далее они будут использованы при построении численной схемы $p-\nu=\text{const}$ (вдоль семейства $C^{+}$ ), $p+\nu=\text{const}$ (вдоль семейства $C^{-}$ ), Результатом первого этапа является получение уравнений характеристик и характеристических соотношений (следствий из них): $\begin{aligned}& {\left( \frac{dx}{dt} \right)\!}_{1}=1 & & (x=t+\text{const}); & & p-\nu=\text{const},\\ & {\left( \frac{dx}{dt} \right)\!}_{2}=-1 & & (x=-t+\text{const}); & & p+\nu=\text{const}.\end{aligned}$ (8.74) К системе (8.74) следует добавить начальные условия $\nu(x,0)=\cos x,~ p(x,0)=2x+1$ . Второй этап. Конструирование численной схемы и вычисление значений искомых функций. 1. Для простоты изложения выделим на оси Ох отрезок $[0;l]$ и на нем зададим три точки $0;0,5;1$ (рис.8.19, а). Они соответствуют нулевому временному слою. Значения функций $p(x,t),\, \nu(x,t)$ в этих точках будем обозначать $p_{i}^0,\,\nu_{i}^0,~ i=0,1,2$ , а сами точки $(x_{i},t^n)$ парами индексов $(i,n)$ , где $i$ — номер точки по переменной $x$ , а $n$ — номер временного слоя. Для $n=0$ с учетом начальных условий $\nu(x,0)=\cos x,~ p(x,0)=2x+1$ получаем $\begin{aligned}& \nu_{0}^{0}= \nu(0;0)=1; & & \nu_{1}^{0}= \nu(0,\!5;0)=0,\!8776; & & \nu_{2}^{0}= \nu(1;0)=0,\!5403;\\[2pt] & p_{0}^{0}= p(0;0)=1; & & p_{1}^{0}= \nu(0,\!5;0)=2; & & p_{2}^{0}= p(1;0)=3. \end{aligned}$ 2. Выведем разностные уравнения, определяющие положение точек на следующем временном слое и значения искомых функций в них. Система (8.74) заменяется системой разностных уравнений. Численное решение задачи строится на характеристической сетке, которая формируется одновременно с процессом нахождения решения. При этом обычно рассматривается трехточечный угловой шаблон с узлами "1","2","3", составленный отрезками характеристик (из семейств $C^{+}$ и $C^{-}$ ), пересекающихся в рассчитываемой точке "3" (рис. 8.19,б). Назовем этот шаблон базовым. Предполагается, что положение узлов "1" , "2" и значения решения в узлах известны из начальных условий или рассчитаны на предыдущем шаге. Положение узла "3" и значение решения в нем подлежат определению. Аппроксимируя производную по формулам ${\left(\frac{dx}{dt}\right)}_{1}\cong \frac{x_3-x_1}{t_3-t_1}$ (вдоль характеристики из семейства $C^{+}$ с $\lambda_1=1$ ), ${\left(\frac{dx}{dt}\right)}_{2}\cong \frac{x_3-x_2}{t_3-t_2}$ (вдоль характеристики из семейства $C^{-}$ с $\lambda_2=-1$ ), получаем $x_3-x_1=t_3-t_1,\qquad x_3-x_2=t_2-t_3.$ (8.75) Условие $p-\nu=\text{const}$ вдоль характеристик из семейств $C^{+}$ записывается для лежащих на ней узлов "1" и "3": $p_3-\nu_3=p_1-\nu_1,$ (8.76) а условие $p-\nu=\text{const}$ вдоль характеристики из семейства $C^{-}$ записывается для лежащих на ней узлов "2" и "3": $p_3+\nu_3= p_2+\nu_2.$ (8.77) Система (8.75)–(8.77) является замкнутой относительно искомых значений $x_3,\,t_3,\,p_3,\,\nu_3$ при известных $x_1, t_1, x_2, t_2, p_1, \nu_1, p_2, \nu_2$ . Формулы (8.75)–(8.77) являются базовыми для расчета всех точек, принадлежащих области $D$ . После совмещения узлов "1" и "2" базового шаблона с выбранными очередными узлами сетки в (8.75)–(8.77) подставляются конкретные величины в соответствии с принятыми обозначениями. 3. Производится расчет первого временного слоя. В решаемой задаче на нулевом временном слое угловой шаблон может занимать два положения (рис. 8.20, а,б). Для точки с номерами $(0,1)$ первого рассчитываемого слоя узлы "1" и "2" углового шаблона совмещаются с узлами $(0,0)$ и $(1,0)$ (рис. 8,20,а). Поэтому система (8.75)–(8.77) имеет вид $\begin{aligned} & x_{0,1}-x_{0,0}= t_{0,1}-t_{0,0}, & & x_{0,1}-x_{1,0}= t_{1,0}-t_{0,1},\\[2pt] & p_0^1-\nu_0^1= p_0^0-\nu_0^0, & & p_0^1+ \nu_0^1= p_1^0+ \nu_1^0. \end{aligned}$ Подставляя $t_{0,0}=0,~ t_{1,0}=0,~ x_{0,0}=0,~ x_{1,0}=0,\!5$ и известные значения $p_0^0, p_1^0, \nu_0^0, \nu_1^0$ на нулевом слое (см. п.1), получаем $\begin{gathered}x_{0,1}-0=t_{0,1}-0;\qquad x_{0,1}-0,\!5=0-t_{0,1};\\[2pt] p_0^1-\nu_0^1=1-1=0;\qquad p_0^1+\nu_0^1=2+ 0,\!8776= 2,\!8776. \end{gathered}$ Решение этой системы: $x_{0,1}=0,\!25;~ t_{0,1}=0,\!25;~ p_0^1=1,\!4388;~ \nu_0^1=1,\!4388$ . Для второй точки первого слоя (рис. 8.20,б) с номерами $(l;l)$ узлы "1" и "2" углового шаблона совмещаются с узлами $(l;0)$ и $(2;0)$ . Тогда система (8.75)–(8.77) имеет вид: $\begin{gathered}x_{1,1}-x_{1,0}=t_{1,1}-t_{1,0};\qquad x_{1,1}-x_{2,0}= t_{2,0}-t_{1,1};\\[2pt] p_1^1-\nu_1^1= p_1^0-\nu_1^0;\qquad p_1^1+ \nu_1^1= p_2^0+ \nu_2^0. \end{gathered}$ Подставляя $t_{1,0}=0,~ t_{2,0}=0,~ x_{1,0}=0,\!5,~ x_{2,0}=1$ и известные значения $p_1^0,\nu_1^0,o_2^0,\nu_2^0$ на нулевом слое (см. п.1), получаем $\begin{gathered}x_{1,1}-0,\!5=t_{11}-0;\qquad x_{1,1}-1=0-t_{1,1};\\[2pt] p_1^1-\nu_1^1=2-0,\!8776=1,\!1224;\qquad p_1^1+\nu_1^1=3+0,\!5403= 3,\!5403. \end{gathered}$ Отсюда $x_{1,1}= 0,\!75;~ t_{1,1}= 0,\!25; p_1^1= 2,\!33135;~ \nu_1^1= 1,\!20895$ . 4. Для расчета точки $(0,2)$ следующего временного слоя узлы "1" и "2" углового шаблона совмещаются с узлами $(0,1)$ и $(1,1)$ (рис.8.20,в). Тогда система (8.75)–(8.77) принимает вид: $\begin{gathered}x_{0,2}-x_{0,1}=t_{0,2}-t_{0,1};\qquad x_{0,2}-x_{1,1}=t_{1,1}-t_{0,2};\\[2pt] p_0^2-\nu_0^2= p_0^1-\nu_0^1;\qquad p_0^2+\nu_0^2= p_1^1+ \nu_1^1. \end{gathered}$ Подставляя значения с предыдущего слоя, имеем $\begin{gathered}x_{0,2}-0,\!25=t_{0,2}-0,\!25;\qquad x_{0,2}-0,\!75= 0,\!25-t_{0,2};\\[2pt] p_0^2-\nu_0^2= 1,\!4388-1,\!4388;\qquad p_0^2+ \nu_0^2=2,\!33135+ 1,\!20895. \end{gathered}$ Отсюда $x_{0,2}= 0,\!5;~ t_{0,2}= 0,\!5;~ p_0^2=\nu_0^2= 1,\!77015$ . Очевидно (см. рис. 8.19,а), количество узлов на каждом последующем слое уменьшается на единицу. Расчет продолжается до замыкания области, так как на последнем слое останется один узел. В данном примере этим узлом является точка $(0,2)$ . Замечание. Рассмотренная в примере 8.10 задача является всего лишь модельной, выбранной для пояснения основных принципов конструирования метода характеристик и его численной реализации. В вычислительной практике, как правило, ставятся нелинейные задачи. Проиллюстрируем некоторые особенности их решения на следующем примере. Пример 8.11. Получить соотношения метода характеристик для решения нелинейной задачи Коши $\begin{gathered}\begin{aligned}&a_{11}\frac{\partial u_1(x,t)}{\partial t}+ a_{12}\frac{\partial u_2(x,t)}{\partial t}+ b_{11}\frac{\partial u_1(x,t)}{\partial x}+ b_{12}\frac{\partial u_2(x,t)}{\partial x}=f_1,\\ &a_{21}\frac{\partial u_1(x,t)}{\partial t}+ a_{22}\frac{\partial u_2(x,t)}{\partial t}+ b_{21}\frac{\partial u_1(x,t)}{\partial x}+ b_{22}\frac{\partial u_2(x,t)}{\partial x}=f_2,\end{aligned}\quad (x,t)\in D\\[4pt] u_1(x,0)= \varphi_1(x),\quad u_2(x,0)= \varphi_2(x),\quad x\in \mathbb{R},\end{gathered}$ где $a_{ij}=a_{ij}(x,t,u),~ b_{ij}=b_{ij}(x,t,u),~ f_{i}=f_{i}(x,t,u),~ i=1,2;~ j=1,2;~ \varphi_{1}(x),~ \varphi_2(x)$ — заданные функции; $D=\mathbb{R}\times (0,T),~ T>0$ — заданный момент времени. [reshenie=Решение, формулы (8.78)-(8.80)] Поставленная задача является частным случаем задачи (8.65) при $m=2$ . Применим метод характеристик, включающий два этапа. Первый этап. Нахождение характеристик и характеристических соотношений. Уравнение (8.67) имеет вид $\|B-A \lambda\|= \begin{vmatrix}b_{11}-\lambda a_{11}& b_{12}-\lambda a_{12}\\ b_{21}-\lambda a_{21}& b_{22}-\lambda a_{22}\end{vmatrix}=0$ или $a\,\lambda^2+ b\,\lambda+c=0,$ где $a=a_{11}a_{12}-a_{12}a_{21},~ b=a_{21}b_{12}+ a_{12}b_{21}-a_{11}b_{22}-a_{22}b_{11},~ c=b_{11}b_{22}-b_{12}b_{21}$ . Решая полученное уравнение, находим выражения для двух характеристических направлений в каждой точке расчетной области: $\begin{aligned}&{\left(\frac{dx}{dt}\right)}_{1}= \lambda_{1}(x,t,u_1,u_2)= \frac{-b+ \sqrt{b^2-ac}}{2a}\quad (C^{+});\\[2pt] &{\left(\frac{dx}{dt}\right)}_{2}= \lambda_{2}(x,t,u_1,u_2)= \frac{-b-\sqrt{b^2-ac}}{2a}\quad (C^{-}). \end{aligned}$ (8.78) Сформируем расширенную матрицу $\bigl(B\,dt-A\,dx,\, f\,dt-A\,du\bigr)= \begin{pmatrix}b_{11}dt-a_{11}dx& b_{12}dt-a_{12}dx& f_1dt-a_{11}du_1-a_{12}du_2\\ b_{21}dt-a_{21}dx& b_{22}dt-a_{22}dx& f_2dt-a_{21}du_1-a_{22}du_2 \end{pmatrix}$ Характеристические соотношения получаются в результате приравнивания нулю определителя второго порядка $(m=2)$ , составленного из второго и третьего столбцов расширенной матрицы: $\begin{vmatrix}b_{12}dt-a_{12}dx& f_1dt-a_{11}du_1-a_{12}du_2\\ b_{22}dt-a_{22}dx& f_2dt-a_{21}du_1-a_{22}du_2 \end{vmatrix}=0.$ Раскрыв определитель и поделив на $dt\ne0$ , получим $\left(e-a\,\frac{dx}{dt}\right)\!du_1+ g\,du_2+ h\,dt+ q\,dx=0,$ где $e=a_{11}b_{22}-a_{21}b_{12},~ g=a_{12}b_{22}-a_{22}b_{12},~ h=b_{12}f_2-b_{22}f_1,~ q=a_{22}f_1-a_{12}f_2$ . Последние два слагаемых перепишем в эквивалентной форме $h\,dt+ q\,dx= \left(h+ q\,\frac{dx}{dt}\right)\!dt$ . Поскольку для семейства характеристик $C^{+}$ выполняется равенство ${\left(\frac{dx}{dt}\right)}_{1}= \lambda_{1}$ , а для семейства $C^{-}$ — равенство ${\left(\frac{dx}{dt}\right)}_{2}= \lambda_{2}$ , то характеристические \соотношения имеют вид $\begin{aligned}& (e-a \lambda_1)\,du_1+ g\,du_2+ (h+q \lambda_1)\,dt=0\quad (C^{+});\\[2pt] & (e-a \lambda_2)\,du_1+ g\,du_2+ (h+q \lambda_2)\,dt=0\quad (C^{-}). \end{aligned}$ (8.79) В результате получена характеристическая система (8.78), (8.79) из четырех обыкновенных дифференциальных уравнений. Они позволяют найти значения $t,x,u_1,u_2$ вдоль характеристик. Второй этап. Конструирование численной схемы. Как и в предыдущем примере, будем использовать трехточечный угловой шаблон (рис. 8.19, б). Предполагается, что в точках "1" и "2" известны значения $t,x,u_1,u_2$ из начальных данных или из расчета предыдущего слоя, т.е. $t_1,x_1,u_{11},u_{21}$ для точки "1" и $t_2,x_2,u_{12},u_{22}$ для точки "2". Здесь при обозначении значений функций $u_1,u_2$ , а далее и при обозначении других функций в точках шаблона первый нижний индекс соответствует номеру самой функции, а второй индекс — номеру точки в шаблоне. Требуется найти значения $t_3,x_3, u_{13}, u_{23}$ определяющие положение точки "3" и значения искомых функций в ней. Предполагается, что точки "1" и "2" близки друг к другу, поэтому заменим дифференциалы конечными приращениями, а производную $\frac{dx}{dt}$ — конечноразностным отношением: $\begin{gathered}dt\cong \Delta t=t_3-t_1,\quad dx\cong \Delta x=x_3-x_1,\quad du_1\cong \Delta u_1=u_{13}-u_{11},\quad du_2\cong \Delta u_2=u_{23}-u_{21},\quad {\left(\frac{dx}{dt} \right)}_{1}\cong \frac{x_3-x_1}{t_3-t_1}\quad (C^{+});\\[4pt] dt\cong \Delta t=t_3-t_2,\quad dx\cong \Delta x=x_3-x_2,\quad du_1\cong \Delta u_1=u_{13}-u_{12},\quad du_2\cong \Delta u_2=u_{23}-u_{22},\quad {\left(\frac{dx}{dt}\right)}_{2}\cong \frac{x_3-x_2}{t_3-t_2}\quad (C^{+}). \end{gathered}$ В итоге уравнения характеристик являются уравнениями прямых, соединяющих точки "1" и "2" с точкой "3", а характеристическая система (8.78), (8.79) преобразуется к виду $\begin{gathered}\begin{aligned}& x_3-x_1= (\lambda_1)_1(t_3-t_1),\\ & x_3-x_2= (\lambda_2)_2(t_3-t_2),\end{aligned}\\[4pt] \begin{aligned}& (e-a \lambda_1)_1(u_{13}-u_{11})+ (g)_1(u_{23}-u_{21})+ (h+ q \lambda_1)_1 (t_3-t_1)=0,\\ & (e-a \lambda_2)_2(u_{13}-u_{12})+ (g)_2(u_{23}-u_{22})+ (h+ q \lambda_2)_2 (t_3-t_2)=0. \end{aligned}\end{gathered}$ (8.80) Здесь нижние индексы 1 и 2 указывают, что значения всех функций, входящих в соответствующие выражения, подсчитываются в точке "1" или "2". Решая систему (8.80), можно найти значения $t_3,x_3,u_{13},u_{23}$ " определяющие положение точки "3" и значения искомых функций в ней. Сначала обычно решаются первые два уравнения и определяются значения $t_3,x_3$ , а затем — последние два уравнения, откуда находятся $u_{13},u_{23}$ . Опишем теперь общую схему численного решения поставленной задачи Коши. 1. При $t=0$ на оси абсцисс задаются точки с некоторым шагом. В этих точках (на рис. 8.21 выделен фрагмент, содержащий точки 1-6) значения функций $u_{1},u_{2}$ определяются начальными данными, т.е. функциями $\varphi_1(x), \varphi_2(x)$ . 2. В соответствии с изложенной методикой последовательно рассчитываются координаты точек 7-11 и приближенные значения решения в них. При этом точки "1" и "2" углового шаблона последовательно совмещаются с точками 1 и 2, затем с точками 2 и 3, 3 и 4, 4 и 5, 5 и 6. Узлами выстраиваемой характеристической сетки служат точки пересечения характеристик (они образуют следующий временной слой). 3. Расчет слоев продолжается до достижения последней точки (точки 21) треугольной области. Отметим, что дифференциальным уравнениям (8.78) соответствуют криволинейные характеристики, а точка "3" должна быть точкой пересечения этих криволинейных характеристик. При переходе от дифференциальных уравнений (8.78),(8.79) к разностным дифференциалы были заменены конечными разностями, а производные конечно-разностными отношениями. Вследствие этого криволинейные характеристики были аппроксимированы отрезками прямых. Решение системы (8.80) объявляется первым приближением, т.е. $t_3^{(1)}= t_3,~ x_3^{(1)}= x_3,$ $u_{13}^{(1)}= u_{13},~ u_{23}^{(1)}= u_{23},$ и может возникнуть необходимость в уточнении координат точки "3" и значений искомых функций в этой точке. Ограничимся изложением одного из возможных способов уточнения. После нахождения первого приближения следует найти средние арифметические значения: $\begin{aligned}& \overline{\lambda}_{1}= \frac{(\lambda_1)_1+ (\lambda_1)_3}{2}\,,\quad \overline{\lambda}_{2}= \frac{(\lambda_2)_2+ (\lambda_2)_3}{2}\,,\quad K_1= \frac{(e-a \lambda_1)_1+ (e-a \lambda_1)_3}{2}\,,\quad K_2= \frac{(e-a \lambda_2)_2+ (e-a \lambda_2)_3}{2}\,,\\[4pt] & \overline{L}_{1}= \frac{(g)_1+ (g)_3}{2}\,,\quad \overline{L}_{2}= \frac{(g)_2+ (g)_3}{2}\,,\quad M_1= \frac{(h+ q \lambda_1)_1+ (h+ q \lambda_1)_3}{2}\,,\quad M_2= \frac{(h+ q \lambda_2)_2+ (h+ q \lambda_2)_3}{2}\,,\end{aligned}$ где нижний индекс 3 указывает, что соответствующее выражение подсчитывается в найденной точке "3". Эти соотношения получаются, если при интегрировании дифференциальных уравнений применять правило трапеций. Затем решается система $\begin{gathered}x_{3}^{(2)}-x_{1}= \overline{\lambda}_{1}(t_{3}^{(2)}-t_{1}),\\ x_{3}^{(2)}-x_{2}= \overline{\lambda}_{2}(t_{3}^{(2)}-t_{2}),\\[4pt] K_1\,\bigl(u_{13}^{(2)}-u_{11}\bigr)+ L_1\,\bigl(u_{23}^{(2)}-u_{21}\bigr)+ M_1\,\bigl(t_{3}^{(2)}-t_{1}\bigr)=0,\\ K_2\,\bigl(u_{13}^{(2)}-u_{12}\bigr)+ L_2\,\bigl(u_{23}^{(2)}-u_{22}\bigr)+ M_2\,\bigl(t_{3}^{(2)}-t_{2}\bigr)=0 \end{gathered}$ относительно нового (второго) приближения $t_2^{(2)},\,x_3^{(2)},\, u_{13}^{(2)},\, u_{23}^{(2)}$ . Описанный процесс продолжают до тех пор, пока решения, полученные при двух последовательных приближениях, не будут совпадать с заданной точностью. Как правило, если расстояние между точками "1" и "2" невелико, достаточно сделать два уточнения. Замечания 1. Приведенные примеры иллюстрируют важное свойство метода характеристик, отличающее его от ранее рассмотренных разностных методов. Оно заключается в том, что расчеты осуществляются вдоль характеристик — линий, по которым распространяется информация об исходных данных, что придает методу четкий физический смысл. 2. Описанный метод расчета, когда точки следующего слоя заранее неизвестны, а вычисления ведутся вдоль характеристик семейств $C^{+}$ и $C^{-}$ , выделяемых в процессе расчета, называется классическим методом характеристик. Существуют модификации этого метода, когда точки следующего слоя задаются заранее, а для нахождения значений искомых функций в них характеристики выпускаются в направлении к известному слою. Такие методы называются сеточно-характеристическими. Метод Годунова В механике сплошных сред при решении задач, описывающихся квазилинейными дифференциальными уравнениями в частных производных, в расчетных областях могут возникать поверхности сильного разрыва, к которым относятся ударные волны, центрированные волны разрежения, поверхности контактного разрыва. На этих поверхностях терпят разрыв как искомые функции, характеризующие исследуемые процессы, так и их производные. При решении таких задач используются два класса численных схем, в одном из которых с использованием законов сохранения явно выделяются поверхности разрыва, а в другом эти поверхности явно не выделяются, а рассчитываются "насквозь" вместе с областями непрерывности. К последнему классу методов относится метод С. К. Годунова. Идею его конструирования рассмотрим на примере решения смешанной краевой задачи для уравнений акустики. Пример 8.12. Получить соотношения метода Годунова для решения задачи $\frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0$ с начальными условиями $u(x,0)=\psi_1(x),~ p(x,0)=\psi_2(x),~ 0 \leqslant x \leqslant 1$ , и краевыми условиями $p(0,t)=\varphi_1(t),~ p(1,t)=\varphi_2(t),~ 0<t \leqslant T$ . Решение, формулы (8.81)-(8.90) Система двух дифференциальных уравнений описывает распространение плоских звуковых волн (малых возмущений) в покоящейся среде. Здесь $u(x,t),\, p(x,t)$ — малые отклонения скорости и давления от их значений в невозмущенной среде, вызванные распространением звуковых волн; $\rho_0$ — плотность невозмущенной среды; $a_0$ — скорость звука. Рассмотрим основные этапы конструирования метода Годунова. 1. Введем разбиение отрезка $[0;1]$ с шагом $h=1\!\!\not{\phantom{\|}}\,n$ , где $n$ — число промежутков. В результате получим узлы $x_{i},~ i=0,1,\ldots,n$ , определяемые соотношением $x_{i+1}=x_{i}+h,\quad i=0,1,\ldots,n-1;\quad x_0=0,\quad x_n=1.$ 2. Заменим функции $\psi_1(x),\, \psi_2(x)$ на начальном слое $t=0$ кусочно-постоянными (пример описанной замены для некоторой произвольной функции $\psi_1(x)$ изображен на рис. 8.22) путем их интегрального осреднения по формулам $u_{1-1\!\not{\phantom{\|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} \psi_1(x)\,dx\,,\quad p_{1-1\!\not{\phantom{\|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} \psi_2(x)\,dx\,,\quad i=1,2,\ldots,n.$ (8.81) Здесь и далее все заданные или ранее рассчитанные функции снабжаются нижними полуцелыми индексами. Из рис. 8.22 видно, что полученная кусочно-постоянная функция имеет в узлах разрывы первого рода. 3. Выберем шаг $\tau$ по времени $t$ и образуем следующий расчетный слой $(t_1=\tau)$ с теми же узлами $x_{i}$ (рис. 8.23). Таким образом, образуется совокупность прямоугольных ячеек. Предполагается, что искомые функции на горизонтальных и вертикальных границах образованных ячеек также являются постоянными. На рассчитываемом слое при $t=\tau$ приближенное решение задачи также представляется в виде кусочно-постоянной функции, имеющей разрывы в узлах сетки, совпадающих с узлами на предыдущем слое. Это означает, что во времени на слоях $t=\text{const}$ сохраняется структура решения, изображенная на рис. 8.22. Значения искомых функций на рассчитываемом слое (верхних границах ячеек) обозначаются верхними полуцелыми индексами: $u^{1-1\!\not{\phantom{\|}}\,\,2},\, p^{1-1\!\not{\phantom{\|}}\,\,2},$ $i=1,2,\ldots,n$ . На вертикальных границах ячеек значения функций обозначаются прописными буквами с нижними целыми индексами: $U_{i},\,P_{i},~ i=0,\ldots,n$ . При $i=0,~ i=n$ значения $P_0$ и $P_n$ определяются краевыми условиями по формулам, аналогичным (8.81): $P_0= \frac{1}{\tau} \int\limits_{0}^{\tau} \varphi_1(t)\,dt\,,\quad P_n= \frac{1}{\tau} \int\limits_{0}^{\tau} \varphi_2(t)\,dt\,.$ (8.82) Значения $U_1,\ldots,U_{n-1},\, P_1,\ldots,P_{n-1}$ находятся в результате решения вспомогательной задачи о распаде произвольного разрыва, описанного в п.4, а значения $U_0,\,U_n$ определяются в п.5. 4. Выберем две соседние ячейки с узлами $x_{i-1},\,x_{i}$ и $x_{i},\,x_{i+1}$ (эти ячейки выделены на рис. 8.23 жирными линиями). С целью нахождения значений $U_{i},\,P_{i}$ на внутренней границе ячеек при $x=x_{i}$ рассмотрим решение вспомогательной задачи, называемой задачей о распаде произвольного разрыва. Дана система уравнений $\frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0$ (8.83) с начальными условиями вида $u(x,0)= \left\{\begin{aligned}&u_{i-1\!\not{\phantom{\|}}\,\,2},\quad x_{i-1} \leqslant x<x_{i},\\ &u_{i+1\!\not{\phantom{\|}}\,\,2},\quad x_{i} \leqslant x<x_{i+1}, \end{aligned}\right.\qquad p(x,0)= \left\{\begin{aligned}&p_{i-1\!\not{\phantom{\|}}\,\,2},\quad x_{i-1} \leqslant x<x_{i},\\ &p_{i+1\!\not{\phantom{\|}}\,\,2},\quad x_{i} \leqslant x<x_{i+1}. \end{aligned}\right.$ Требуется найти решение $u(x,t),\,p(x,t)$ при $t>0$ и $x=x_{i}$ , то есть $u(x_{i},t)=U_{i},~ p(x_{i},t)=P_{i}$ . Заметим, что в данной задаче начальные данные имеют разрыв первого рода в точке $x_{i}$ . При $t>0$ они трансформируются, т.е. происходит распад разрыва, что обусловливает название рассматриваемой задачи. Для ее решения применим метод характеристик. Перепишем систему (8.83) в матричной форме $A\cdot \frac{\partial u(x,t)}{\partial t}+ B\cdot \frac{\partial u(x,t)}{\partial x}=0$ , где $A=\begin{pmatrix}1&0\\ 0&1 \end{pmatrix}\!,~~ B=\begin{pmatrix}0& 1\!\!\not{\phantom{\|}}\,\rho_0\\ \rho_0a_0^2&0 \end{pmatrix}\!,~~ u=\begin{pmatrix}u\\p \end{pmatrix}$ . Для нахождения характеристических направлений используем условие $\|B-A\cdot \lambda\|= \begin{vmatrix}-\lambda& 1\!\!\not{\phantom{\|}}\,\rho_0\\ \rho_0a_0^2&-\lambda\end{vmatrix}=0$ . Отсюда ${\left(\frac{dx}{dt}\right)\!}^{+}= \lambda^{+}= a_0,~~ {\left(\frac{dx}{dt}\right)\!}^{-}= \lambda^{-}=-a_0$ . Интегрируя последние уравнения, имеем $x-a_0t=\text{const}~(C^{+}),\qquad x+a_0t=\text{const}~(C^{-}).$ Если при $t=0$ выполняется условие $x=x_{i}$ , то из последних соотношений выделяются две прямые, исходящие из точки $x_{i}$ на начальном слое: $x-a_0t=x_{i}~(C^{+}),\qquad x+a_0t=x_{i}~(C^{-}).$ Для получения характеристических соотношений запишем расширенную матрицу (8.69) (учитывая, что $f\equiv 0$ ): $\left(B\,dt-A\,dx,-A\! \begin{pmatrix}du\\dp\end{pmatrix}\right)= \begin{pmatrix}-dx&\dfrac{1}{\rho_0}\,dt&-dt\\ \rho_0a_0^2dt&-dx&-dp \end{pmatrix}\!.$ Приравнивая определитель матрицы, составленной из двух последних столбцов, нулю, получаем $-\dfrac{1}{\rho_0}\,dp\,dt-dx\,du=0.$ Полагая $dt\ne0$ , разделим уравнение на $dt\colon\,-\dfrac{1}{\rho_0}\,dp-\frac{dx}{dt}\,du=0.$ Подставляя вместо $\frac{dx}{dt}$ сначала $\lambda^{+}=a_0$ , а затем $\lambda^{-}=-a_0$ , получаем $du+ \dfrac{1}{\rho_0\,a_0}\,dp=0\quad (C^{+}),\qquad du-\dfrac{1}{\rho_0\,a_0}\,dp= 0\quad (C^{-}).$ Так как коэффициенты обоих уравнений постоянны, то их можно переписать в эквивалентной форме $d\! \left(u+ \dfrac{1}{\rho_0\,a_0}\,p\right)=0\quad (C^{+}),\qquad d\! \left(u-\dfrac{1}{\rho_0\,a_0}\,p\right)=0\quad (C^{-}).$ Интегрируя и обозначая постоянные интегрирования $I^{+},\,I^{-}$ , имеем $u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}\quad (C^{+}),\qquad u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}\quad (C^{-}).$ (8.84) Полученные соотношения называются инвариантами Римана. Они показывают, что величина $u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}$ остается постоянной вдоль семейства $C^{+}$ , а величина $u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}$ остается постоянной вдоль семейства $C^{-}$ . Теперь, используя направления характеристик и сохраняющиеся вдоль них инварианты, построим численную схему решения задачи о распаде разрыва. Рассмотрим две характеристики $x-a_0t=x_{i-1}$ и $x-a_0t=x_{i}$ из семейства $C^{+}$ , исходящие из узлов $x_{i-1},\,x_{i}$ , а также две характеристики $x+a_0t=x_{i}$ и $x+a_0t=x_{i+1}$ из семейства $C^{-}$ , исходящие из узлов $x_{i},\,x_{i+1}$ (рис. 8.24). В результате пересечения указанных характеристик образуются области I, II, III , первые две из которых представлены треугольниками, а третья — четырехугольником. Вертикальная граница $(x_{i}= \text{const})$ , разделяющая ячейки, на которой ищется решение задачи, располагается в области III. Покажем, что на этой границе искомое решение постоянно. Для этого выберем на ней произвольную точку (назовем ее точка 3) , а решение в ней обозначим $U_{i},\, P_{i}$ . В эту точку приходит характеристика семейства $C^{+}$ , исходящая из точки 1 промежутка $(x_{i-1}, x_{i})$ . Точке 1 соответствуют заданные значения $u_{i-1\!\not{\phantom{\|}}\,\,2},\, p_{i-1\!\not{\phantom{\|}}\,\,2}$ . Вдоль этой характеристики остается постоянной величина $u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}$ , то есть $U_{i}+ \frac{1}{\rho_0\,a_0}\cdot P_{i}= u_{i-1\!\not{\phantom{\|}}\,\,2}+ \frac{1}{\rho_0\,a_0}\cdot p_{i-1\!\not{\phantom{\|}}\,\,2}.$ Кроме того, в точку 3 приходит характеристика семейства $C^{-}$ , исходящая из точки 2 промежутка $(x_{i},x_{i+1})$ . Точке 2 соответствуют заданные значения и $u_{i+ 1\!\not{\phantom{\|}}\,\,2},\, p_{i+ 1\!\not{\phantom{\|}}\,\,2}$ . Вдоль этой характеристики остается постоянной величина $u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}$ , то есть $U_{i}-\frac{1}{\rho_0\,a_0}\cdot P_{i}= u_{i+ 1\!\not{\phantom{\|}}\,\,2}-\frac{1}{\rho_0\,a_0}\cdot p_{i+ 1\!\not{\phantom{\|}}\,\,2}.$ Разрешая систему из двух последних уравнений относительно $U_{i}$ и $P_{i}$ , получаем $\begin{gathered}U_{i}= \frac{u_{i-1\!\not{\phantom{\|}}\,\,2}+ u_{i+ 1\!\not{\phantom{\|}}\,\,2}}{2}-\frac{p_{i+ 1\!\not{\phantom{\|}}\,\,2}-p_{i+ 1\!\not{\phantom{\|}}\,\,2}}{2\rho_0\,a_0}\,\\ P_{i}= \frac{p_{i-1\!\not{\phantom{\|}}\,\,2}+ p_{i+ 1\!\not{\phantom{\|}}\,\,2}}{2}+ \rho_0\,a_0\cdot \frac{u_{i+ 1\!\not{\phantom{\|}}\,\,2}-u_{i-1\!\not{\phantom{\|}}\,\,2}}{2}\,,\end{gathered}\quad i=1,2,\ldots,n-1.$ (8.85) Поскольку точка 3 была выбрана произвольно на вертикальной границе ячейки, то полученные выражения определяют решение в любой точке этой границы (очевидно, решение постоянно вдоль границы), за исключением самой точки $x_{i}$ и точки пересечения характеристик, исходящих из точек $x_{i-1},\, x_{i+1}$ . Таким образом, формулы (8.85) определяют "большие" величины $U_{i},\, P_{i}$ на всех внутренних вертикальных границах ячеек. Результаты решения данной вспомогательной задачи определяют шаг $\tau$ по времени. Действительно, влияние каждого из разрывов, относящихся к соседним узлам $x_{i-1},\,x_{i}$ на слое $t=0$ , ограничивается характеристиками, исходящими из этих узлов. Физически это означает, что разрывы, возникшие при $t=0$ , распространяются со скоростью звука $a_0$ . Характеристики, исходящие из соседних узлов навстречу друг другу, пересекутся в точке с координатами $t=\frac{h}{2a_0},~ x=\frac{x_{i-1}+ x_{i}}{2}$ . От этой встречи разрывов возникнут новые волны, которые через следующий отрезок времени длительностью $\frac{h}{2a_0}$ достигнут вертикальных границ $(x_{i-1}=\text{const},~ x_{i}=\text{const})$ соседних ячеек. Именно до этого момента, т.е. в течение времени $\frac{h}{a_0}= \frac{h}{2a_0}+ \frac{h}{2a_0}$ , значения $U_{i},\, P_{i}$ будут сохраняться неизменными и определяться по формулам (8.85). Таким образом, величина шага по времени должна выбираться из условия $\tau< \frac{h}{a_0}$ . При $i=0$ для нахождения величины $U_0$ при известном $P_0$ используется инвариант Римана $u-\frac{p}{\rho_0 a_0}=I^{-}$ вдоль семейства $C^{-}$ (рис. 8.25). Поскольку на промежутке $(0,x_1)$ заданы значения $u_{1\!\not{\phantom{\|}}\,\,2},\,p_{1\!\not{\phantom{\|}}\,\,2}$ , то $U_0-\frac{P_0}{\rho_0a_0}= u_{1\!\not{\phantom{\|}}\,\,2}-\frac{p_{1\!\not{\phantom{\|}}\,\,2}}{\rho_0 a_0}$ . Отсюда получаем $U_0= \frac{1}{\rho_0a_0}\cdot P_0+ u_{1\!\not{\phantom{\|}}\,\,2}-\frac{1}{\rho_0 a_0}\cdot p_{1\!\not{\phantom{\|}}\,\,2}.$ (8.86) При $i=n$ для нахождения величины $U_n$ при известном $P_n$ используется инвариант Римана $u+\frac{p}{\rho_0a_0}=I^{+}$ вдоль семейства $C^{+}$ (рис. 8.25). Поскольку на последнем промежутке $(x_{n-1},1)$ заданы значения $u_{n-1\!\not{\phantom{\|}}\,\,2},\, p_{n-1\!\not{\phantom{\|}}\,\,2}$ , то $U_n+\frac{P_n}{\rho_0a_0}= u_{n-1\!\not{\phantom{\|}}\,\,2}+ \frac{P_{n-1\!\not{\phantom{\|}}\,\,2}}{\rho_0a_0}$ . Отсюда получаем $U_n=-\frac{1}{\rho_0a_0}\cdot P_n+ u_{n-1\!\not{\phantom{\|}}\,\,2}+ \frac{1}{\rho_0a_0}\cdot p_{n-1\!\not{\phantom{\|}}\,\,2}.$ (8.87) На рис. 8.25 показаны две крайние ячейки с тремя характеристиками, одна из которых (верхняя) ограничивает область влияния начальных данных, заданных на промежутках $(x_0,x_1)$ и $(x_{n-1},x_{n})$ соответственно. 6. Перейдем к задаче расчета искомых величин на следующем временном слое. Для этого рассмотрим одну ячейку (рис. 8.26), для которой известны значения приближенного решения $u_{i-1\!\not{\phantom{\|}}\,\,2},\, p_{i-1\!\not{\phantom{\|}}\,\,2}$ на нижней границе и $U_{i-1},P_{i-1},U_{i},P_{i},$ $i=1,2,\ldots,n$ , на левой и правой вертикальных границах (см. п.4). Требуется определить значения $p^{i-1\!\not{\phantom{\|}}\,\,2},\, u^{i-1\!\not{\phantom{\|}}\,\,2}$ на верхней границе ячейки. Для решения поставленной задачи перепишем уравнения (8.83) в так называемом дивергентном виде $\frac{\partial u(x,t)}{\partial t}+ \frac{\partial}{\partial x}\! \left(\frac{p(x,t)}{\rho_0} \right)=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \frac{\partial}{\partial x}\bigl(\rho_0a_0^2u(x,t) \bigr)=0.$ (8.88) Запишем двойной интеграл от левой и правой частей уравнений по поверхности ячейки. Для его нахождения применим формулу Грина перехода от двойного интеграла к криволинейному интегралу второго рода: $\iint\limits_{S}\! \left(\frac{\partial \Phi_1}{\partial t}+ \frac{\partial \Phi_2}{\partial x}\right)\!dx\,dt= \oint\limits_{\gamma}\bigl(\Phi_1\,dx-\Phi_2\,dt\bigr),$ где $S$ — поверхность ячейки, а $\gamma$ — ее граница, обходимая против часовой стрелки. Таким образом, вместо системы (8.88) получается система интегральных уравнений $\oint\limits_{\gamma}u\,dx-\frac{p}{\rho_0}\,dt=0,\qquad \oint\limits_{\gamma}p\,dx-\rho_0a_0^2u\,dt=0.$ Проводя интегрирование, имеем $\begin{gathered}u_{i-1\!\not{\phantom{\|}}\,\,2}h-\frac{P_{i}}{\rho_0}\cdot \tau-u^{i-1\!\not{\phantom{\|}}\,\,2}h+ \frac{P_{i-1}}{\rho_0}\cdot \tau=0,\\[4pt] p_{i-1\!\not{\phantom{\|}}\,\,2}h-\rho_0a_0^2U_{i}\tau-p^{i-1\!\not{\phantom{\|}}\,\,2}h+ \rho_0a_0^2U_{i-1}\tau=0. \end{gathered}$ Из этих уравнений получаем $\begin{gathered}u^{i-1\!\not{\phantom{\|}}\,\,2}= u_{i-1\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h\,\rho_0}\cdot (P_{i}-P_{i-1}),\\[4pt] p^{i-1\!\not{\phantom{\|}}\,\,2}= p_{i-1\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h}\cdot \rho_0a_0^2(U_{i}-U_{i-1}). \end{gathered}$ Таким образом, получены формулы для расчета первого временного слоя $(t=\tau)$ . После нахождения значений $u^{\frac{1}{2}},p^{\frac{1}{2}}, u^{\frac{3}{2}}, p^{\frac{3}{2}},\ldots, u^{n-\frac{1}{2}}, p^{n-\frac{1}{2}}$ из условия $\tau<\frac{h}{a_0}$ выбирается очередной шаг и формируется новый временной слой $(t_2=t_1+ \tau)$ . Найденные значения используются в качестве опорных (они переобозначаются путем присвоения им соответствующих нижних индексов) для нахождения приближенного решения на следующем слое по формулам (8.85)–(8.87), (8.89). Заметим, что для вычисления $P_0,\,P_n$ вместо (8.82) используются формулы интегрального осреднения, аналогичные (8.81). Для расчета, например, второго временного слоя они имеют вид $P_0= \frac{1}{\tau} \int\limits_{t_1}^{t_2} \varphi_1(t)\,dt\,,\qquad P_n= \frac{1}{\tau} \int\limits_{t_1}^{t_2} \varphi_2(t)\,dt\,.$ (8.90) Изложенная схема метода Годунова имеет первый порядок аппроксимации по $\tau$ и $h$ . Замечание. Из процедуры интегрирования следует, что искомые значения $u^{n-1\!\not{\phantom{\|}}\,\,2},\, p^{n-1\!\not{\phantom{\|}}\,\,2}$ соответствуют среднеинтегральным величинам $u^{n-1\!\not{\phantom{\|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} u(x,\tau)\,dx\,,\qquad p^{n-1\!\not{\phantom{\|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} p(x,\tau)\,dx\,.$ Приведем пример расчета по описанной схеме для частного случая рассматриваемой задачи (где $\rho_0=0,\!123;~ a_0=0,\!25$ ): $\begin{gathered}\begin{aligned}& \frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\\ & \frac{\partial p(x,t)}{\partial t}+\rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0,\end{aligned}\qquad 0<x<1,\quad 0<t<1,\\[4pt] \begin{aligned}& u(x,0)= \psi_1(x)=\sin x, & & p(x,0)=\psi_2(x)= \rho_0x\sin x, & & 0 \leqslant x \leqslant 1,\\ & p(0,t)= \varphi_1(t)=0,\!1t+ \rho_0a_0^2\sin t, & & p(1,t)= \varphi_2(t)= 0,\!1t+ \rho_0a_0^2\sin t, & & 0<t \leqslant1, \end{aligned} \end{gathered}$ Решение А. Выберем число промежутков $n=2$ . Тогда находится шаг $h=\frac{1}{2}= 0,\!5$ и вводятся узлы $x_0=0;~ x_1=0,\!5;~ x_2=1$ . Зададим шаг $\tau= 0,\!5$ по времени из условия $\tau< \frac{h}{a_0}= \frac{0,\!5}{0,\!25}=2$ . Сформируем две ячейки "A" и "B" первого временного слоя $t_1=\tau= 0,\!5$ (рис. 8.27). Б. Заменим функции $\psi_1(x)=\sin x,~ \psi_2(x)=\rho_0x\sin x$ кусочно-постоянными по формулам (8.81). Для вычисления интегралов используются квадратурные формулы трапеций: $u_{i-1\!\not{\phantom{\|}}\,\,2}= \frac{1}{h}\cdot \frac{u_{i-1}+ u_{i}}{2}\cdot h=\frac{u_{i-1}+ u_{i}}{2}\,,\qquad p_{i-1\!\not{\phantom{\|}}\,\,2}= \frac{1}{h}\cdot \frac{p_{i-1}+ p_{i}}{2}\cdot h=\frac{p_{i-1}+ p_{i}}{2}\,,\quad i=1;2,$ где $u_{i}=\sin x_{i},~ p_{i}=\rho_0x_{i}\sin x_{i}$ . Поэтому $\begin{aligned}& u_{1\!\not{\phantom{\|}}\,\,2}\cong \frac{0+0,\!479}{2}\cong 0,\!239; & & u_{3\!\not{\phantom{\|}}\,\,2}\cong \frac{0,\!479+ 0,\!841}{2}\cong 0,\!660;\\[2pt] <br />& p_{1\!\not{\phantom{\|}}\,\,2}\cong \frac{0+0,\!0294}{2}\cong 0,\!239; & & p_{3\!\not{\phantom{\|}}\,\,2}\cong \frac{0,\!479+ 0,\!841}{2}\cong 0,\!660<br />\end{aligned}$ В. Найдем значения $P_0,\,P_2$ по формулам (8.82). Для вычисления интегралов, как и ранее, применяется формула трапеций: $P_0\cong \frac{0+0,\!0537}{2}\cong 0,\!0269;\qquad P_2\cong \frac{0+0,\!0463}{2}\cong 0,\!0231.$ Рассчитаем значение $U_0$ на левой границе ячейки "А" по формуле (8.86): $U_0=\frac{P_0}{\rho_0a_0}+ u_{1\!\not{\phantom{\|}}\,\,2}-\frac{p_{1\!\not{\phantom{\|}}\,\,2} }{\rho_0a_0}\cong 0,\!635$ и значение $U_2$ на правой границе ячейки по формуле (8.87): $U_2=-\frac{P_2}{\rho_0a_0}+ u_{3\!\not{\phantom{\|}}\,\,2}+ \frac{p_{3\!\not{\phantom{\|}}\,\,2} }{\rho_0a_0}\cong 1,\!476$ . Г. Вычислим значения $U_1,\,P_1$ на правой границе ячейки "А" (левой границе ячейки "В" ). По формулам (8.85) получаем $U_1=\frac{u_{1\!\not{\phantom{\|}}\,\,2}+ u_{3\!\not{\phantom{\|}}\,\,2}}{2}-\frac{p_{3\!\not{\phantom{\|}}\,\,2}-p_{1\!\not{\phantom{\|}}\,\,2}}{2\rho_0a_0}\cong-0,\!392;\qquad \frac{p_{1\!\not{\phantom{\|}}\,\,2}+ p_{3\!\not{\phantom{\|}}\,\,2}}{2}-\rho_0a_0\, \frac{u_{3\!\not{\phantom{\|}}\,\,2}-u_{1\!\not{\phantom{\|}}\,\,2}}{2}\cong 0,\!0341.$ Д. Вычислим значения $u^{1\!\not{\phantom{\|}}\,\,2},\, p^{1\!\not{\phantom{\|}}\,\,2},\, u^{3\!\not{\phantom{\|}}\,\,2},\, p^{3\!\not{\phantom{\|}}\,\,2}$ на верхних границах ячеек: $\begin{aligned}& u^{1\!\not{\phantom{\|}}\,\,2}= u_{1\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_1-P_0)\cong 0,\!1805; & & p^{1\!\not{\phantom{\|}}\,\,2}= p_{1\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h}\,\rho_0a_0^2(U_1-U_0)\cong 0,\!0226;\\[4pt] & u^{3\!\not{\phantom{\|}}\,\,2}= u_{3\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_2-P_1)\cong 0,\!749; & & p^{3\!\not{\phantom{\|}}\,\,2}= p_{3\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h}\,\rho_0a_0^2(U_2-U_1)\cong 0,\!0521. \end{aligned}$ Таким образом, первый слой полностью рассчитан. Е. Для расчета второго слоя $(t_2=t_1+\tau)$ с шагом $\tau=0,\!5$ полученные значения $u^{1\!\not{\phantom{\|}}\,\,2},\, p^{1\!\not{\phantom{\|}}\,\,2},\, u^{3\!\not{\phantom{\|}}\,\,2},\, p^{3\!\not{\phantom{\|}}\,\,2}$ переобозначим, для чего верхние индексы заменяются нижними: $u_{1\!\not{\phantom{\|}}\,\,2}=0,\!1805;\quad p_{1\!\not{\phantom{\|}}\,\,2}= 0,\!0226;\quad u_{3\!\not{\phantom{\|}}\,\,2}=0,\!749;\quad p_{3\!\not{\phantom{\|}}\,\,2}= 0,\!0521$ Сформируем две ячейки "С" и "D" второго временного слоя (см. рис. 8.27). Ж. Найдем значения $P_0,\,P_2$ по формулам (8.82). Для вычисления интегралов, как и ранее, применяется формула трапеций: $P_0\cong \frac{0,0537+ 0,1065}{2}\cong 0,\!0801;\qquad P_2\cong \frac{0,0463+ 0,00935}{2}\cong 0,\!0278.$ Рассчитаем значение $U_0$ на левой границе ячейки "С" по формуле (8.86): $U_0= \frac{1}{\rho_0\,a_0}\cdot P_0+ u_{1\!\not{\phantom{\|}}\,\,2}-\frac{1}{\rho_0\, a_0}\cdot p_{1\!\not{\phantom{\|}}\,\,2}\cong 2,\!0504$ и значение $U_2$ на правой границе ячейки "D" по формуле (8.87): $U_2=-\frac{1}{\rho_0\,a_0}\cdot P_2+ u_{3\!\not{\phantom{\|}}\,\,2}+ \frac{1}{\rho_0\,a_0}\cdot p_{3\!\not{\phantom{\|}}\,\,2}\cong 0,\!8651.$ 3. Вычислим значения $U_1,\,P_1$ на правой границе ячейки "С" (левой границе ячейки "D"). По формулам (8.85) получаем $U_1= \frac{u_{1\!\not{\phantom{\|}}\,\,2}+ u_{3\!\not{\phantom{\|}}\,\,2}}{2}-\frac{p_{3\!\not{\phantom{\|}}\,\,2}-p_{1\!\not{\phantom{\|}}\,\,2}}{2\,\rho_0\,a_0}\cong-0,\!0150;\qquad P_1=\frac{p_{1\!\not{\phantom{\|}}\,\,2}+ p_{3\!\not{\phantom{\|}}\,\,2}}{2}-\rho_0a_0\, \frac{u_{3\!\not{\phantom{\|}}\,\,2}-u_{1\!\not{\phantom{\|}}\,\,2}}{2}\cong 0,\!02861.$ Вычислим значения $u^{1\!\not{\phantom{\|}}\,\,2},\, p^{1\!\not{\phantom{\|}}\,\,2},\, u^{3\!\not{\phantom{\|}}\,\,2},\, p^{3\!\not{\phantom{\|}}\,\,2}$ на верхних границах ячеек "С" и "D": $\begin{aligned}& u^{1\!\not{\phantom{\|}}\,\,2}= u_{1\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_1-P_0)\cong 0,\!599; & & p^{1\!\not{\phantom{\|}}\,\,2}= p_{1\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h}\,\rho_0\,a_0^2(U_1-U_0)\cong 0,\!0385;\\[2pt] & u^{3\!\not{\phantom{\|}}\,\,2}= u_{3\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_2-P_1)\cong 0,\!7556; & & p^{3\!\not{\phantom{\|}}\,\,2}= p_{3\!\not{\phantom{\|}}\,\,2}-\frac{\tau}{h}\,\rho_0\,a_0^2(U_2-U_1)\cong 0,\!04533. \end{aligned}$ Таким образом, второй слой полностью рассчитан и найдено приближенное решение поставленной задачи. Замечание. Здесь рассмотрен пример методики конструирования метода С.К. Годунова и ее реализации для относительно простой линейной задачи. В нелинейном случае метод существенно усложняется. Для более глубокого изучения проблемы авторы рекомендуют, где данный метод применяется для различных классов нелинейных задач и задач с тремя и четырьмя независимыми переменными. Математический форум (помощь с решением задач, обсуждение вопросов по математике). Если заметили ошибку, опечатку или есть предложения, напишите в комментариях.

Численные методы решения уравнений в частных производных

Онлайн-сервисы

Нахождение НОД и НОК Разложение числа на простые множители Сравнения по модулю Операции над множествами Операции над векторами Разложение вектора по базису. Доказательство, что векторы образуют базис Чертёж треугольника по координатам вершин Решение треугольника Решение Пирамиды Построение Пирамиды по координатам вершин Чертёж многоугольника по координатам вершин Решение систем методом Крамера и Матричным Онлайн построение графика кривой 2-го порядка Определение вида кривой или поверхности 2-го порядка по инвариантам МНК и регрессионный анализ Онлайн + графики Онлайн число, сумма и дата прописью

Алгоритмы JavaScript

Алгоритмы поиска Алгоритмы сортировки Уникальные элементы массива Объединение, пересечение и разность массивов НОД и НОК Операции над матрицами Дата прописью

Введение в анализ

Функции: понятие, определение, графики Непрерывность функции Исследование функции и построение графика

Теория множеств

Множества: понятие, определение, примеры Точечные множества Замкнутые и открытые множества Мера множества Группы, кольца, поля в математике Поле комплексных чисел Кольцо многочленов Основная теорема алгебры и ее следствия

Математическая логика

Алгебра высказываний

Аксиоматика и логические рассуждения Методы доказательств теорем Алгебра высказываний и операции над ними Формулы алгебры высказываний Тавтологии алгебры высказываний Логическая равносильность формул Нормальные формы для формул высказываний Логическое следование формул Приложение алгебры высказываний для теорем Дедуктивные и индуктивные умозаключения Решение логических задач Принцип полной дизъюнкции

Булевы функции

Множества, отношения и функции в логике Булевы функции от одного и двух аргументов Булевы функции от n аргументов Системы булевых функций Применение булевых функций к релейно-контактным схемам Релейно-контактные схемы в ЭВМ Практическое применение булевых функций

Теория формального

Формализованное исчисление высказываний Полнота и другие свойства формализованного исчисления высказываний Независимость системы аксиом формализованного исчисления высказываний

Логика предикатов

Логика предикатов Логические операции над предикатами Кванторные операции над предикатами Формулы логики предикатов Тавтологии логики предикатов Преобразования формул и следование их предикатов Проблемы разрешения для общезначимости и выполнимости формул Применение логики предикатов в математике Строение математических теорем Аристотелева силлогистика и методы рассуждений Принцип полной дизъюнкции в предикатной форме Метод полной математической индукции Необходимые и достаточные условия Логика предикатов и алгебра множеств Формализованное исчисление предикатов

Неформальные и формаль-
ные аксиоматические теории

Неформальные аксиоматические теории Свойства аксиоматических теорий Формальные аксиоматические теории Формализация теории аристотелевых силлогизмов Свойства формализованного исчисления предикатов Формальные теории первого порядка Формализация математической теории

Теория алгоритмов

Интуитивное представление об алгоритмах Машины Тьюринга и тезис Рекурсивные функции Нормальные алгоритмы Маркова Разрешимость и перечислимость множеств Неразрешимые алгоритмические проблемы Теорема Гёделя о неполноте формальной арифметики

Математическая логика и компьютеры

Математическая логика и языки программирования Применение компьютеров для доказательства теорем математической логики Математическая логика и логическое программирование Математическая логика и информатика Математическая логика и искусственный интеллект

Дискретная математика

Множества и отношения

Теория множеств: понятия и определения Операции над множествами Кортеж и декартово произведение множеств Соответствия и бинарные отношения на множествах Операции над соответствиями на множествах Семейства множеств Специальные свойства бинарных отношений Отношения эквивалентности на множестве Упорядоченные множества Теорема о неподвижной точке Мощность множества Парадокс Рассела Метод характеристических функций

Группы и кольца

Алгебраические структуры и операции Группоиды, полугруппы, группы Кольца, тела, поля Области целостности в теории колец Модули и линейные пространства Подгруппы и подкольца Теорема Лагранжа о порядке конечной группы Гомоморфизмы групп и нормальные делители Гомоморфизмы и изоморфизмы колец Алгебра кватернионов

Полукольца и булевы алгебры

Полукольца: определение, аксиомы, примеры Замкнутые полукольца Полукольца и системы линейных уравнений Булевы алгебры и полукольца Решетки и полурешетки

Алгебраические системы

Алгебраические системы: модели и алгебры Подсистемы алгебраических систем Конгруэнции и фактор-системы Гомоморфизмы алгебраических систем Прямые произведения алгебраических систем Конечные булевы алгебры Многосортные алгебры

Теория графов

Теория графов: основные понятия и определения Способы представления графов Неориентированные и ориентированные деревья Остовное дерево и алгоритм Краскала Методы систематического обхода вершин графа Алгоритмы поиска в глубину и ширину в графах Задача о путях во взвешенных ориентированных графах Изоморфизм, гомоморфизм и автоморфизм графов Топологическая сортировка вершин графа Элементы цикломатики в теории графов

Булева алгебра и функции

Булевы функции и булев куб Таблицы булевых функций и булев оператор Равенство булевых функций. Фиктивные переменные Формулы и суперпозиции булевых функций Дизъюнктивные и конъюнктивные нормальные формы Построение минимальных ДНФ Теорема Поста и классы Критерий Поста Схемы из функциональных элементов

Конечные автоматы и регулярные языки

Конечные автоматы и регулярные языки Алфавит, слово, язык в программировании Порождающие грамматики (грамматики Хомского) Классификация грамматик и языков Регулярные языки и регулярные выражения Конечные автоматы Допустимость языка конечным автоматом Теорема Клини Детерминизация конечных автоматов Минимизация конечных автоматов Лемма о разрастании для регулярных языков Обоснование алгоритма детерминизации автоматов Конечные автоматы с выходом Морфизмы и конечные подстановки Машины Тьюринга

Контекстно-свободные языки

Контекстно-свободные языки и грамматики Приведенная форма КС-грамматики Лемма о разрастании для КС-языков Магазинные автоматы (автомат с магазинной памятью) Алгоритм построения МП-автомата по КС-грамматике Алгоритм построения КС-грамматики по МП-автомату Алгебраические свойства КС-языков Основное свойство суперпозиции КС-языков Пересечение контекстно-свободных языков Методы синтаксического анализа КС-языков Восходящий синтаксический анализ и LR(k)-грамматики Семантика формальных языков Принцип индукции по неподвижной точке Графовое представление МП-автоматов

Интегральное исчисление

Неопределённый и определённый

Неопределенный и определенный интегралы Свойства интегралов Интегрирование по частям Интегрирование методом замены переменной Интегрирование различных рациональных функций Интегрирование различных иррациональных функций Интегрирование различных тригонометрических функций Определенный интеграл и его основные свойства Необходимое и достаточное условие интегрируемости Теоремы существования первообразной Свойства определенных интегралов Несобственные интегралы Интегральное определение логарифмической функции

Приложения интегралов

Вычисление площадей плоских фигур Площади фигур в различных координатах Вычисление объемов тел с помощью интегралов Объём тела вращения Вычисление длин дуг кривых Формулы длины дуги регулярной кривой Кривизна плоской кривой Площадь поверхности вращения тела

Интегралы в физике

Статические моменты и координаты центра тяжести Теоремы Гульдина–Паппа Вычисление моментов инерции Другие приложения интегралов в физике

Основные интегралы

Интеграл Ньютона-Лейбница Интеграл Римана Интеграл Лебега

Вариационное исчисление

Примеры вариационных задач Дифференциальное уравнение Эйлера Функционалы, зависящие от нескольких функций Задача о минимуме кратного интеграла

Финансовый анализ

Анализ эффективности

Критерии и показатели эффективности предприятия Методы анализа эффективности деятельности Факторный анализ прибыли от операционной деятельности Анализ безубыточности предприятия Операционный рычаг и эффект финансового рычага Анализ и оценка состава, структуры и динамики доходов и расходов Анализ рентабельности и резервов устойчивого роста капитала Анализ распределения прибыли предприятия Анализ и оценка чувствительности показателей эффективности

Анализ устойчивости

Финансовая устойчивость и долгосрочная платежеспособность Характеристика типов финансовой устойчивости

Рыночная активность

Финансовый анализ рыночной активности Методика анализа рыночной активности Анализ и оценка дивидендного дохода на одну акцию

Инвестиционная деятельность

Инвестиции: экономическая сущность и классификация Государственное регулирование инвестиционной деятельности Источники финансовых ресурсов на капитальные вложения Инвестиции в основные фонды Оценка состояния основных фондов Амортизация основных фондов Капитальное строительство в инвестиционном процессе Планирование инвестиций в форме капитальных вложений Экономическая эффективность инвестиций Финансирование капитальных вложений Кредитование капитальных вложений Кредитоспособность Финансирование и кредитование затрат Финансирование и кредитование инвестиционной деятельности потребительской кооперации Финансирование и кредитование капитальных вложений потребительской кооперации Инвестиционное строительное проектирование

Анализ инвестиций

Инвестиции и инвестиционная деятельность предприятия Задачи финансового анализа инвестиций предприятия Учет фактора времени в инвестиционной деятельности Аннуитет и финансовая рента в инвестициях Учет фактора инфляции при инвестировании Оценка фактора риска инвестиционного проекта Методы оценки эффективности инвестиций Показатели эффективности инвестиционного проекта

Стоимость компании

Концепция построения международных стандартов финансовой отчетности (МСФО) Экономическое содержание международных стандартов финансовой отчётности Цели и принципы оценки стоимости акций и активов компании Оценка акций и активов предприятия по справедливой стоимости Методы оценки справедливой стоимости акций предприятия Затратный подход к оценки стоимости компаний и акций Сравнительный подход к оценки стоимости предприятий и акций Доходный подход к оценке стоимости компании и акций Выбор ставки дисконтирования при инвестировании в акции Метод капитализации прибыли Сравнение подходов к оценке стоимости компаний и пакетов акций

Форвардные контракты

Форвардный контракт и цена Форвардная цена акции на бирже Цена форвардного контракта инвестора Форвардная цена акции с учетом величины дивиденда Форвардная цена акции с учетом ставки дивиденда Форвардная цена валюты на рынке форекс Форвардный валютный курс и инфляция на рынке Форвардная цена товара и спотовый рынок Форвардная цена при различии ставок по кредитам и депозитам Синтетический форвардный контракт на акции и валюту

Теория вероятностей

Основные понятия теории вероятностей Зависимые и независимые случайные события Повторные независимые испытания Формула Бернулли Одномерные случайные величины Многомерные случайные величины Функции случайных величин Законы распределения целочисленных случайных величин Законы распределения непрерывных случайных величин Предельные теоремы теории вероятностей Закон больших чисел и предельные теоремы Вероятностные закономерности

Математическая статистика

Элементы математической статистики Выборочный метод Оценки параметров генеральной совокупности Статистические гипотезы Критерии согласия Теоретические и эмпирические частоты

Теория очередей (СМО)

Определение системы массового обслуживания Уравнения Колмогорова Предельные вероятности состояний Определение СМО с отказами Определение СМО с ожиданием (очередью)

Аналитическая геометрия

Векторная алгебра

Метрические понятия и аксиомы геометрии Равенство и подобие геометрических фигур Бинарные отношения Вектор, его направление и длина Линейные операции над векторами Линейная зависимость и независимость векторов Отношение коллинеарных векторов Проекции векторов на прямую и на плоскость Угол между векторами Ортогональные проекции векторов Координата вектора на прямой и базис Координаты вектора на плоскости и базис Координаты вектора в пространстве и базис Операции над векторами в координатной форме Ортогональный и ортонормированный базисы Cкалярное произведение векторов и его свойства Выражение скалярного произведения через координаты векторов Векторное произведение векторов и его свойства Смешанное произведение векторов и его свойства Ориентированные площади и объемы Двойное векторное произведение и его свойства Применение векторов в задачах на аффинные свойства фигур Применение произведений векторов при решении геометрических задач Применение векторной алгебры в механике

Системы координат

Прямоугольные координаты Преобразования прямоугольных координат Полярная система координат Цилиндрическая система координат Сферические координаты Аффинные координаты Аффинные преобразования координат Аффинные преобразования плоскости Примеры аффинных преобразований плоскости Аффинные преобразования пространства Многомерное координатное пространство Линейные и аффинные подпространства Скалярное произведение n-мерных векторов Преобразования систем координат

Геометрия на плоскости

Алгебраические линии на плоскости Общие уравнения геометрических мест точек Алгебраические уравнения линий на плоскости Уравнения прямой, проходящей через точку перпендикулярно вектору Уравнения прямой, проходящей через точку коллинеарно вектору Уравнения прямой, проходящей через две точки Уравнения прямой с угловым коэффициентом Взаимное расположение прямых Примеры задач с прямыми на плоскости Системы неравенств с двумя неизвестными Системы линейных уравнений с двумя неизвестными

Линии 2-го порядка

Канонические уравнения линий второго порядка Порядок приведения уравнения линии к каноническому виду Эллипс Гипербола Парабола Квадратичные неравенства с двумя неизвестными Применение линий 1-го и 2-го порядков в задачах на экстремум функций

Инварианты линий

Классификация линий 2-го порядка по инвариантам Приведение уравнения линии к каноническому виду по инвариантам

Геометрия в пространстве

Способы задания ГМТ в пространстве Алгебраические уравнения поверхностей Уравнения плоскости, проходящей через точку перпендикулярно вектору Уравнения плоскости, компланарной двум неколлинеарным векторам Уравнения плоскости, проходящей через три точки Взаимное расположение плоскостей Типовые задачи с плоскостями Уравнения прямых в пространстве Взаимное расположение прямых в пространстве Типовые задачи с прямыми в пространстве

Поверхности 2-го порядка

Канонические уравнения поверхностей Порядок приведения уравнения поверхности к каноническому виду Поверхности второго порядка Эллипсоиды Гиперболоиды Конусы Параболоиды Применение поверхностей 1-го и 2-го порядков в задачах на экстремум функций

Инварианты поверхностей

Классификация поверхностей 2-го порядка по инвариантам Квадратичные неравенства с тремя неизвестными Приведение уравнения поверхности к канониче-скому виду по инвариантам

Линейная алгебра

Матрицы и операции

Линейные операции над матрицами Умножение матриц Возведение матриц в степень Многочлены от матриц Транспонирование и сопряжение матриц Блочные матрицы Произведение и сумма матриц Кронекера Метод Гаусса приведения матрицы к ступенчатому виду Элементарные преобразования матриц

Определители

Определители матриц и их основные свойства Формула полного разложения определителя Формула Лапласа полного разложения определителя Определитель произведения матриц Методы вычисления определителей

Ранг матрицы

Линейная зависимость и линейная независимость строк (столбцов) матрицы Ранг матрицы и базисный минор матрицы Методы вычисления ранга матрицы Ранг системы столбцов (строк)

Обратная матрица

Обратные матрицы и их свойства Ортогональные и унитарные матрицы Способы нахождения обратной матрицы Матричные уравнения Односторонние обратные матрицы Скелетное разложение матрицы Полуобратная матрица Псевдообратная матрица

Системы уравнений

Системы линейных алгебраических уравнений Метод Гаусса решения систем линейных уравнений Структура общего решения системы уравнений Решение систем с помощью полуобратных матриц Псевдорешения системы линейных уравнений

Функциональные матрицы

Функциональные матрицы скалярного аргумента Производные матриц по векторному аргументу Линейные и квадратичные формы и их преобразования Приведение форм к каноническому виду Закон инерции вещественных квадратичных форм Знакоопределенность форм вещественных квадратичных Формы и исследование функций на экстремум

Многочленные матрицы

Многочленные матрицы (лямбда-матрицы) Операции над лямбда-матрицами Простые преобразования многочленных матриц Инвариантные множители многочленной матрицы

Функции от матриц

Собственные векторы и значения матрицы Подобие числовых матриц Характеристический многочлен матрицы Минимальный многочлен матрицы Теорема Гамильтона-Кэли Жорданова форма матрицы Приведение матрицы к жордановой форме Многочлены от матриц Применение многочленов от матриц Функции от матриц

Линейные пространства

Линейные пространства: определение и примеры Линейная зависимость и независимость n-мерных векторов Размерность и базис линейного пространства Преобразования координат в линейном пространстве Изоморфизм линейных пространств

Подпространства

Подпространства линейного пространства Пересечение и сумма подпространств Способы описания подпространств Нахождение дополнения и суммы подпространств Нахождение пересечения подпространств

Линейные отображения

Линейные многообразия Линейные отображения Матрица линейного отображения Ядро и образ линейного отображения

Линейные операторы

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

Евклидовы пространства

Комплексный анализ

Комплексные числа

Комплексные числа в алгебраической форме Комплексные числа в тригонометрической и показательной формах Множества на комплексной плоскости Последовательности и ряды комплексных чисел

Комплексные функции

Функции комплексного переменного. Предел, непрерывность и производная Элементарные функции комплексного переменного Дифференцирование функций комплексного переменного Аналитические функции и их свойства Конформные отображения
и их свойства Интегрирование функций комплексного переменного

Функциональные ряды в комплексной области

Функциональные ряды и последовательности Степенные ряды и их свойства Разложение функций в степенные ряды Нули аналитических функций Ряд Лорана и разложение функций по целым степеням

Особые точки, Вычеты

Изолированные особые точки функций и полюсы Вычеты и их применение Вычисление интегралов с помощью вычетов Вычеты и расположение нулей многочлена

Операционное исчисление

Преобразование Лапласа и его свойства Решение ДУ операционным методом Анализ выходных процессов линейных стационарных систем Z-преобразование и его свойства

Дифференциальные уравнения

ДУ первого порядка

Основные понятия и определения ДУ Метод изоклин для ДУ 1-го порядка Метод последовательных приближений ДУ с разделяющимися переменными Однородные ДУ Линейные ДУ 1-го порядка Дифференциальное уравнение Бернулли ДУ в полных дифференциалах Интегрирующий множитель ДУ, не разрешенные относительно производной Дифференциальное уравнение Риккати Составление ДУ семейств линий Задачи на траектории Особые решения ДУ

ДУ высших порядков

Понятия и определения ДУ высших порядков ДУ, допускающие понижение порядка Линейная независимость функций Определители Вронского и Грама Однородные и неоднородные дифференциальные уравнения Задача Коши и Уравнение Эйлера Линейные ДУ с переменными коэффициентами Метод Лагранжа решения ДУ Краевые задачи для ДУ высших порядков Разложение решения ДУ в степенной ряд Разложение решения ДУ в обобщенный степенной ряд Нахождение периодических решений ДУ Асимптотическое интегрирование ДУ

Системы ДУ

Системы ДУ: понятия и определения Сведение системы ДУ к одному уравнению Нахождение интегрируемых комбинаций Интегрирование однородных линейных систем ДУ Методы интегрирования неоднородных систем ДУ Преобразование Лапласа и решение ДУ и систем

Теория устойчивости

Устойчивость решений ДУ по Ляпунову Простейшие типы точек покоя Метод функций Ляпунова Устойчивость решений ДУ по первому приближению Критерии устойчивости Рауса–Гурвица и Михайлова ДУ с малым параметром при производной

Численные методы

Методы алгебры

Численные методы линейной алгебры Численные методы решения СЛАУ Итерационный метод Шульца обратной матрицы Методы решения задач о собственных значениях и векторах матрицы Методы решения нелинейных уравнений Методы решения систем нелинейных уравнений

Методы теории приближений

Методы приближения сеточных функций Методы функциональной интерполяции Методы интегрально-дифференциальной интерполяции Методы интегрального сглаживания Методы интерполяции и сглаживания сплайнами Методы численного дифференцирования и интегрирования Методы численного дифференцирования Методы численного интегрирования

Методы решения обыкновенных ДУ

Численные методы решения задачи Коши Разностные схемы для решения задачи Коши Составные схемы для решения задачи Коши Экстраполяционные методы решения задачи Коши Непрерывно-дискретные методы решения задачи Коши Численные методы решения краевых задач

Методы решения ДУ в частных производных

Численные методы решения уравнений математической физики с двумя переменными Принципы построения разностных схем для уравнений в частных производных Разностные схемы решения уравнений в частных производных 1-го порядка Разностные схемы решения уравнений в частных производных 2-го порядка Численные методы решения уравнений в частных производных Численные методы решения уравнений математической физики с тремя переменными

Численные методы решения уравнений в частных производных

Метод прямых

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

Основную идею конструирования метода прямых рассмотрим на примере нахождениия решения $u(x,y)$ уравнения:

$\begin{gathered}A(x,y)\frac{\partial^2 u(x,y)}{\partial x^2}+ 2B(x,y)\frac{\partial^2 u(x,y)}{\partial y\,\partial x}+ C(x,y) \frac{\partial^2 u(x,y)}{\partial y^2}+ D(x,y)\frac{\partial u(x,y)}{\partial x}\,+\\[2pt] +\,E(x,y)\frac{\partial u(x,y)}{\partial y}+ F(x,y)u(x,y)= G(x,y),\quad (x,y)\in \Omega= (0,L)\times (0,M), \end{gathered}$

(8.59)

с краевыми условиями

$\begin{aligned}& (x,0)= \varphi(x), & & u(x,M)= \psi(x),\\ & (0,y)= \Phi(y), & & u(L,y)= \Psi(y) \end{aligned}$

(8.60)

на границах области прямоугольной формы (рис. 8.15), где $\varphi(x),\psi(x), \Phi(x), \Psi(x)$ — заданные функции. Пусть выполнено условие $B^2-AC<0$ эллиптичности уравнения для всех $(x,y)\in \Omega$ . Отсюда следует, что $A(x,y)\ne 0,~ C(x,y)\ne 0$ .

Нахождение решения ДУ в частных производных на границах области

Для решения задачи выберем на отрезке $[0,M]$ точки $y_{j}= j\,h,~ 0 \leqslant j \leqslant J,~ J\,h=M$ , и проведем через них прямые, параллельные оси $Ox$ . На каждой прямой дифференциальное уравнение в частных производных (8.59) приближенно заменим обыкновенным дифференциальным уравнением для функций $u(x,y_{j})$ . Для этого аппроксимируем частные производные по $y\colon$

$\begin{gathered}\left.{\frac{\partial \widehat{u}(x,y)}{\partial y}}\right|_{y=y_{j}}= \frac{\widehat{u}(x,y_{j+1})-\widehat{u}(x,y_{j-1})}{2h}\,,\\[2pt] \left.{\frac{\partial^2 \widehat{u}(x,y)}{\partial y^2}}\right|_{y=y_{j}}= \frac{\widehat{u}(x,y_{j+1})-2\,\widehat{u}(x,y_{j})+ \widehat{u}(x,y_{j-1})}{h^2}\,,\\[2pt] \left.{\frac{\partial^2 \widehat{u}(x,y)}{\partial y\, \partial x}}\right|_{y=y_{j}}= \frac{\widehat{u}_{x}(x,y_{j+1})-\widehat{u}_{x}(x,y_{j-1})}{2h}\,. \end{gathered}$

(8.61)

Введем обозначения

$\begin{gathered}\widehat{u}(x,y_{j})= \widehat{u}_{j}(x),\qquad \frac{\partial \widehat{u}(x,y_{j})}{\partial x}= \widehat{u}\,'_{j}(x),\qquad \frac{\partial^2 \widehat{u}(x,y_{j}) }{\partial x^2}= \widehat{u}\,''_{j}(x),\\[4pt] A(x,y_{j})= A_{j}(x),\qquad B(x,y_{j})= B_{j}(x),\qquad C(x,y_{j})= C_{j}(x)\quad \text{etc.}\end{gathered}$

(8.62)

На прямых, описываемых уравнениями $y_{j}=\text{const},~ j=1,\ldots,J-1$ , получаем

$\begin{gathered} A_{j}(x)\widehat{u}\,''_{j}(x)+ B_{j}(x)\frac{\widehat{u}\,'_{j+1}(x)-\widehat{u}\,'_{j-1}(x)}{h}+ C_{j}(x) \frac{\widehat{u}_{j+1}(x)-2\,\widehat{u}_{j}(x)+ \widehat{u}_{j-1}(x)}{h^2}+ D_{j}(x) \widehat{u}\,'_{j}(x)+\\[2pt] E_{j}(x) \frac{\widehat{u}_{j+1}(x)-\widehat{u}_{j-1}(x)}{2h}+ F_{j}(x) \widehat{u}_{j}(x)= G_{j}(x),\quad j=1,\ldots,J-1;~ 0<x<L. \end{gathered}$

Здесь $\widehat{u}_{1}(x),\ldots, \widehat{u}_{J-1}(x)$ — приближенное решение задачи на прямых $y=y_1=\text{const},\ldots, y=y_{J-1}= \text{const},$ $A_{j}(x),\ldots, G_{j}(x)$ — коэффициенты уравнения (8.59) на этих прямых. В силу краевых условий находим

$\begin{gathered}\widehat{u}_{0}(x)= u(x,0)= \varphi(x),\qquad \widehat{u}_{J}(x)= u(x,M)= \psi(x),\\[4pt] \widehat{u}_{j}(0)= \Phi(y_{j})= \Phi(j\,h),\quad \widehat{u}_{j}(L)= \Psi(y_{j})= \Psi(j\,h),\quad j=1,\ldots,J-1. \end{gathered}$

(8.63)

Таким образом, от линейного дифференциального уравнения в частных производных (8.59) и краевыми условиями (8.60) переходим к дифференциально-разностной краевой задаче (8.62),(8.63) для системы из $(n-1)$ обыкновенных дифференциальных уравнений (8.62) относительно $(n-1)$ неизвестных функций $\widehat{u}_{1}(x),\ldots, \widehat{u}_{J-1}(x)$ . В результате ее решения можно найти законы изменения искомой функции $u(x,y)$ вдоль прямых $y_{j}=\text{const},$ $j=0,1,\ldots,J$ , то есть $\widehat{u}(x,y_{j}),~ j=0,1,\ldots,J$ . Значения функции в промежуточных точках по у может быть найдено с помощью интерполяции. Система (8.62),(8.63) называется системой уравнений метода прямых.

Замечания

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

2. Метод прямых удобнее применять в том случае, когда коэффициенты в уравнении (8.59) не зависят от $x$ . Тогда система (8.62) будет системой с постоянными коэффициентами.

3. Если расчетная область отлична от прямоугольника, при реализации метода возникают трудности, связанные с тем, что функции $u_{j}(x)$ определены на различных множествах.

4. При решении нестационарных уравнений дискретизация проводится по пространственным переменным, а время $t$ остается непрерывной переменной.

5. На практике система дифференциальных уравнений (8.62) из-за большого числа $J$ должна решаться численно, поэтому потребуется вводить шаг не только по $y$ .

6. Если изменение искомой функции внутри множества $\Omega$ в направлении $y$ не велико, то можно ограничиться одной или двумя прямыми.

Пример 8.9. Составить уравнения метода прямых для решения одномерной задачи теплопроводности (где $\psi(x)$ — заданная функция)

$\begin{aligned}& u_{t}=a^2\cdot u_{xx},\quad 0<x<1,\quad 0<t<1,\\ & u(x,0)= \psi(x),\quad 0 \leqslant x \leqslant 1,\\ & u(0,t)= u(1,t)=0,\quad 0<t \leqslant 1. \end{aligned}$

Решение, формула (8.64)

Для решения задачи выберем на отрезке $[0,1]$ точки $x_{i}= i\,h,~ i=0,1,\ldots, I;~ I\,h=1$ и проведем через них прямые, параллельные оси $Ot$ (рис. 8.16). На каждой прямой приближенно заменим уравнение теплопроводности обыкновенным дифференциальным уравнением для функции $u(x_{i},t)$ . Для этого заменим частную производную $u_{xx}$ по формуле

$\left.{\frac{\partial^2 \widehat{u}}{\partial x^2}}\right|_{x=x_{i}}= \frac{\widehat{u}(x_{i-1},t)-2\,\widehat{u}(x_{i},t)+ \widehat{u}(x_{i+1},t)}{h^2}$

и введем обозначения $\widehat{u}(x_{i},t)= \widehat{u}_{i}(t),~ \frac{\partial \widehat{u}(x_{i}, t)}{\partial t}= \widehat{u}\,'_{i}(t)$ .

Тогда на прямых $x_{i}=\text{const},~ i=1,\ldots,I-1$ , получаем систему и

$\begin{gathered}\widehat{u}\,'_{i}(t)= a^2\cdot \frac{\widehat{u}(x_{i-1},t)-2\,\widehat{u}(x_{i},t)+ \widehat{u}(x_{i+1},t)}{h^2}\,,\quad 1 \leqslant i \leqslant I-1,\\[4pt] \begin{aligned}& \widehat{u}_{0}(t)= u(0,t)=0, & & 0<t \leqslant 1,\\ & \widehat{u}_{I}(t)= u(1,t)=0, & & 0<t \leqslant 1 \end{aligned} \end{gathered}$

(8.64)

с начальными условиями $\widehat{u}_{i}(0)= \psi(x_{i})= \psi(i\,h),~ i=1,2,\ldots,I-1$ . Полученную систему дифференциальных уравнений можно записать в векторно-матричной форме:

$\frac{d}{dt}\! \begin{pmatrix}\widehat{u}_1\\ \vdots\\ \vdots\\ \widehat{u}_{I-1} \end{pmatrix}= \frac{a^2}{h^2}\! \begin{pmatrix}-1&1&\cdots& 0&0\\ 1&-2&1&\cdots&0\\ \vdots& \vdots& \vdots& \vdots& \vdots\\ \vdots& \vdots& 1&-2&1\\ 0&0&\cdots& 1&-2 \end{pmatrix}\!\cdot\! \begin{pmatrix}\widehat{u}_1\\ \vdots\\ \vdots\\ \widehat{u}_{I-1}\end{pmatrix}\!.$

Задача имеет аналитическое решение

$\widehat{u}(t)= H^{-1}\cdot\! \begin{pmatrix}\exp \dfrac{\lambda_1a^2t}{h^2}&0&0\\ 0&\ddots&0\\ 0&0& \exp \dfrac{\lambda_{I-1}a^2t}{h^2} \end{pmatrix}\!\cdot H\,u(0),$

где $u(0)= \begin{pmatrix}\psi(h)\\ \psi(2h)\\ \vdots\\ \psi \bigl((I-1)h\bigr)\end{pmatrix}\!,~ \lambda_{i}=-\left(2\cos \frac{\pi\,i}{n}+ 2\right)\!,~ i=1,\ldots,I-1$ — собственные значения матрицы системы, элементы матрицы $H$ (собственные векторы матрицы $A$ ) равны $H_{ik}= \sin \frac{\pi\,ik}{n}\,,$ $i=1,\ldots,I-1,$ $k=1,\ldots,I-1$ .

Заметим, что при численном решении системы обыкновенных дифференциальных уравнений при увеличении числа $I$ (улучшении точности) коэффициент жесткости $S$ системы увеличивается. При больших значениях $n$ справедливо $S= \frac{\max_{i}|\lambda_{i}|}{\min_{i}|\lambda_{i}|}\cong \frac{4I^2}{\pi^2}$ .

Метод характеристик

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

Рассмотрим основные идеи метода на примере задачи Коши для квазилинейной системы дифференциальных уравнений m-го порядка с двумя независимыми переменными $x,\,t\colon$

$\begin{aligned}&\begin{aligned}&a_{11} \frac{\partial u_{1}(x,t)}{\partial t}+ \ldots+a_{1m} \frac{\partial u_{m}(x,t)}{\partial t}+ b_{11} \frac{\partial u_{1}(x,t)}{\partial x}+\ldots+b_{1m} \frac{\partial u_{m}(x,t)}{\partial x}= f_1(x,t,u),\\[-2pt] &\quad\vdots\\[-4pt] &a_{m1} \frac{\partial u_{1}(x,t)}{\partial t}+ \ldots+a_{mm} \frac{\partial u_{m}(x,t)}{\partial t}+ b_{m1} \frac{\partial u_{1}(x,t)}{\partial x}+\ldots+b_{mm} \frac{\partial u_{m}(x,t)}{\partial x}= f_m(x,t,u),\end{aligned}\\[4pt] &\begin{aligned}& u_1(x,0)= \varphi_1(x), & & x\in \mathbb{R},\\[-2pt] & \quad\vdots\\[-4pt] & u_m(x,0)= \varphi_m(x), & & x\in \mathbb{R}, \end{aligned} \end{aligned}$

где $a_{ij}=a_{ij}(x,t,u),~ b_{ij}=b_{ij}(x,t,u),~ i=1,\ldots,m;~ j=1,\ldots,m;~ D=R\times (0,T),~ T>0$ — заданный момент времени. Перепишем систему в векторно-матричном виде

$\begin{gathered}A(x,t,u)\frac{\partial u(x,t)}{\partial t}+ B(x,t,u)\frac{\partial u_{1}(x,t)}{\partial x}= f(x,t,u),\quad (x,t)\in D,\\[2pt] u(x,0)= \varphi(x),\quad x\in \mathbb{R}, \end{gathered}$

(8.65)

где $A(x,t,u)= \bigl(a_{ij}(x,t,u)\bigr),~ B(x,t,u)= \bigl(b_{ij}(x,t,u)\bigr)$ — матрицы размеров $(m\times m),~ f(x,t,u)= \begin{pmatrix}f_1(x,t,u)\\\vdots\\ f_m(x,t,u) \end{pmatrix}$ — вектор-функция, $u(x,t)= \begin{pmatrix}u_1(x,t)\\ \vdots\\ u_m(x,t) \end{pmatrix}$ — искомая вектор-функция правых частей, $\varphi(x)= \begin{pmatrix} \varphi_1(x)\\ \vdots\\ \varphi_m(x) \end{pmatrix}$ — заданная вектор-функция.

Пусть известно, что система (8.65) имеет гладкое решение в области $D$ и на некоторой кривой $\gamma$ оно задано и равно $\Bigl.{u(x,t)}\Bigr|_{\gamma}$ (рис. 8.17). Рассмотрим задачу о нахождении решения в окрестности кривой $\gamma$ .

Выберем на кривой $\gamma$ произвольную точку $(x_0,t_0)$ , а вектор бесконечно малого смещения вдоль этой кривой из точки $(x_0,t_0)$ обозначим $(dx,dt)$ . Присоединим к уравнениям (8.65) выражение для полного дифференциала (для упрощения записи опустим аргументы):

$A\,\frac{\partial u}{\partial t}+ B\,\frac{\partial u}{\partial x}=f,\qquad du= \frac{\partial u}{\partial x}\,dx+ \frac{\partial u}{\partial t}\,dt\,.$

Умножим первое уравнение на $dt$ , а из второго выразим $\frac{\partial u}{\partial t}\,dt\colon$

$A\,\frac{\partial u}{\partial t}\,dt+ B\,\frac{\partial u}{\partial x}\,dt= f\,dt,\qquad \frac{\partial u}{\partial t}\,dt= du-\frac{\partial u}{\partial x}\,dx\,.$

Подставляя второе уравнение в первое, имеем $A\,du-A\,\frac{\partial u}{\partial x}\,dx+ B\,\frac{\partial u}{\partial x}\,dt= f\,dt$ или

$(B\,dt-A\,dx)\cdot \frac{\partial u}{\partial x}= f\,dt-A\,du\,.$

(8.66)

Если определитель матрицы $(B\,dt-A\,dx)$ отличен от нуля в каждой точке кривой $\gamma$ , то частная производная $\frac{\partial u}{\partial x}$ и, следовательно, частная производная $\frac{\partial u}{\partial t}$ определяются однозначно. В этом случае можно найти решение в окрестности кривой $\gamma\colon\, \Bigl.{u(x,t)= u(x,t)}\Bigr|_{\gamma}+ du$ , что соответствует разложению искомого решения по формуле Тейлора до членов первого порядка. Если определитель матрицы $(B\,dt-A\,dx)$ в каждой точке кривой $\gamma$ равен нулю, то в силу предположения о существовании решения производные $\frac{\partial u}{\partial x},\, \frac{\partial u}{\partial t}$ и, следовательно, функция $u(x,t)$ , находятся неоднозначно. В этом случае кривая у называется характеристикой. Приравнивая определитель матрицы $(B\,dt-A\,dx)$ к нулю и предполагая, что $dt\ne0$ вдоль кривой $\gamma$ , получаем соотношение

$\left|B-A\cdot \frac{\partial u}{\partial t}\right|=0\quad \text{ili}\quad |B-A \lambda|=0,$

(8.67)

где $\lambda(x,t,u)= \frac{\partial u}{\partial t}$ . Линии, задаваемые дифференциалами смещения $dx,\,dt$ , вдоль которых справедливо (8.67), являются характеристиками. Анализ разрешимости уравнения (8.67) позволяет ввести классификацию систем вида (8.65).

Если уравнение (8.67) имеет $m$ различных действительных корней, то система (8.65) называется гиперболической; если их число меньше $m$ , то система называется параболической; если уравнение не имеет действительных корней, то — эллиптической. В дальнейшем будем рассматривать только гиперболические системы, т.е. имеющие $m$ различных характеристик:

$\left(\frac{dx}{dt}\right)_{i}= \lambda_{i}(x,t,u),\quad i=1,2,\ldots,m,$

(8.68)

где функции $\lambda_{i}(x,t,u)$ задают характеристические направления в каждой точке кривой $\gamma$ .

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

$\bigl(B\,dt-A\,dx,\, f\,dt-A\,du\bigr).$

(8.69)

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

Этот способ получения характеристических соотношений основан на применении теоремы Кронекера-Капелли: решение системы (8.66) существует тогда и только тогда, когда ранг матрицы $(B\,dt-A\,dx)$ этой системы равен рангу расширенной матрицы. Предполагается, что решение существует и определитель матрицы системы (m-го порядка) вдоль характеристик равен нулю. Поэтому должен быть равен нулю определитель расширенной матрицы, составленный из произвольных $m$ столбцов.

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

Пример 8.10. Получить соотношения метода характеристик для решения линейной задачи Коши

$\frac{\partial^2 \varphi(x,t)}{\partial t^2} = \frac{\partial^2 \varphi(x,t)}{\partial x^2}\,,\quad (x,t)\in \mathbb{R}\times (0;T),\quad \varphi(x,0)=x^2+x,\quad \varphi_t(x,0)=\cos x,\quad x\in \mathbb{R}.$

Решение, формулы (8.70)-(8.77)

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

Для применения метода характеристик сначала сведем задачу к виду (8.65), т.е. к системе двух уравнений первого порядка. Для этого введем новые искомые функции

$\nu(x,t)= \frac{\partial \varphi(x,t)}{\partial t}\,\qquad p(x,t)= \frac{\partial \varphi(x,t)}{\partial x}\,.$

Здесь $\nu(x,t),\,p(x,t)$ — отклонения скорости и давления от их значений в невозмущенной среде, вызванные распространением звуковых волн, $\varphi(x,t)$ — потенциал функции $\nu(x,t)$ . Тогда уравнение можно представить в форме

$\frac{\partial }{\partial t}\! \left(\frac{\partial \varphi(x,t)}{\partial t}\right)= \frac{\partial }{\partial x}\! \left(\frac{\partial \varphi(x,t)}{\partial x}\right)\quad \text{ili}\quad \frac{\partial \nu(x,t)}{\partial t}-\frac{\partial p(x,t)}{\partial x}=0.$

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

$\frac{\partial^2 \varphi(x,t)}{\partial x\,\partial t}= \frac{\partial^2 \varphi(x,t)}{\partial t\,\partial x}\quad \text{ili}\quad \frac{\partial p(x,t)}{\partial t}= \frac{\partial \nu(x,t)}{\partial x}\,.$

В результате от исходной задачи для волнового уравнения переходим к задаче для системы двух дифференциальных уравнений с начальными условиями относительно новых искомых функций:

$\left\{\begin{aligned}& \frac{\partial \nu(x,t)}{\partial t}-\frac{\partial p(x,t)}{\partial x}=0, & & \nu(x,0)=\cos x,\\ & \frac{\partial p(x,t)}{\partial t}-\frac{\partial \nu(x,t)}{\partial x}=0, & & p(x,0)=2x+1. \end{aligned}\right.$

Запишем систему в виде (8.65)

$A\cdot \frac{\partial u}{\partial t}+ B\cdot \frac{\partial u}{\partial x}=0,\quad u(x,0)= \begin{pmatrix}\cos x\\ 2x+1\end{pmatrix}\!,$

(8.70)

где $A=\begin{pmatrix}1&0\\0&1\end{pmatrix}\!,~ B=\begin{pmatrix}0&-1\\-1&0\end{pmatrix}\!,~ u=\begin{pmatrix}\nu\\p \end{pmatrix}\!,~ m=2$ . Для решения задачи (8.70) применим метод характеристик, содержащий два этапа.

Первый этап. Нахождение характеристик и вывод характеристических соотношений. Для этого решим уравнение (8.67), записанное для (8.70):

$|B-A\cdot \lambda|= \begin{vmatrix}-\lambda&-1\\-1&-\lambda \end{vmatrix}=0.$

Получаем $\lambda^2-1=0,~ \lambda_{1,2}=\pm1$ , т.е. два направления характеристик:

${\left(\frac{dx}{dt}\right)\!}_{1}= \lambda_1=1;\quad {\left(\frac{dx}{dt} \right)\!}_{2}= \lambda_2=-1.$

(8.71)

Интегрируя, найдем общие решения полученных дифференциальных уравнений. Для первого уравнения $x=t+\text{const}$ , а для второго $x=-t+\text{const}$ . Им соответствуют два семейства характеристик (в данной задаче прямых линий). Их принято обозначать буквой $C^{+}$ с верхним индексом "+" для семейства характеристик, имеющих меньший угол с положительным направлением оси Ох по сравнению с соответствующим углом для характеристик семейства $C^{-}$ . В данной задаче характеристики семейства $C^{+}$ имеют положительный наклон, а характеристики семейства $C^{-}$ — отрицательный наклон (рис. 8.18). Условно их также называют характеристиками первого и второго семейств соответственно.

Для получения характеристических соотношений запишем расширенную матрицу (в решаемой задаче $f=\begin{pmatrix}0\\0 \end{pmatrix}$ ):

$\bigl(B\,dt-A\,dx,\, f\,dt-A,\,du\bigr)= \left(\underbrace{\begin{pmatrix}0&-dt\\-dt&0 \end{pmatrix}}_{B\,dt}-\underbrace{\begin{pmatrix}dx&0\\ 0&dx \end{pmatrix}}_{A\,dx},-\underbrace{\begin{pmatrix}1&0\\ 0&1 \end{pmatrix}\!\! \begin{pmatrix}d\nu\\ dp\end{pmatrix}}_{A\,du} \right)= \begin{pmatrix}-dx&-dt-d\nu\\-dt&-dx&-dp \end{pmatrix}\!.$

Приравнивая определитель матрицы, составленной из последних двух столбцов, нулю, получаем

$\begin{vmatrix}-dt&-d\nu\\-dx&-dp \end{vmatrix}= dp\cdot dt-d\nu\cdot dx=0.$

Поделив на $dt$ (предполагается, что $dt\ne 0$ ), имеем $dp-d\nu\, \frac{dx}{dt}=0$ . Так как на прямых семейства $C^{+}$ выполняется равенство ${\left( \frac{dx}{dt} \right)\!}_{1}= \lambda_1=1$ , то соответствующие характеристическое соотношение имеет вид

$dp-d\nu=0.$

(8.72)

Аналогично на прямых семейства $C^{-}$ справедливо ${\left( \frac{dx}{dt} \right)\!}_{2}= \lambda_2=-1$ и

$dp+ d\nu=0.$

(8.73)

Обыкновенные дифференциальные уравнения (8.71)–(8.73) образуют систему уравнений метода характеристик. Полученные характеристические соотношения (8.72), (8.73) в данной задаче могут быть проинтегрированы. В результате получаются конечные соотношения, называемые инвариантами Римана. Далее они будут использованы при построении численной схемы

$p-\nu=\text{const}$ (вдоль семейства $C^{+}$ ),

$p+\nu=\text{const}$ (вдоль семейства $C^{-}$ ),

Результатом первого этапа является получение уравнений характеристик и характеристических соотношений (следствий из них):

$\begin{aligned}& {\left( \frac{dx}{dt} \right)\!}_{1}=1 & & (x=t+\text{const}); & & p-\nu=\text{const},\\ & {\left( \frac{dx}{dt} \right)\!}_{2}=-1 & & (x=-t+\text{const}); & & p+\nu=\text{const}.\end{aligned}$

(8.74)

К системе (8.74) следует добавить начальные условия $\nu(x,0)=\cos x,~ p(x,0)=2x+1$ .

Второй этап. Конструирование численной схемы и вычисление значений искомых функций.

1. Для простоты изложения выделим на оси Ох отрезок $[0;l]$ и на нем зададим три точки $0;0,5;1$ (рис.8.19, а). Они соответствуют нулевому временному слою. Значения функций $p(x,t),\, \nu(x,t)$ в этих точках будем обозначать $p_{i}^0,\,\nu_{i}^0,~ i=0,1,2$ , а сами точки $(x_{i},t^n)$ парами индексов $(i,n)$ , где $i$ — номер точки по переменной $x$ , а $n$ — номер временного слоя. Для $n=0$ с учетом начальных условий $\nu(x,0)=\cos x,~ p(x,0)=2x+1$ получаем

$\begin{aligned}& \nu_{0}^{0}= \nu(0;0)=1; & & \nu_{1}^{0}= \nu(0,\!5;0)=0,\!8776; & & \nu_{2}^{0}= \nu(1;0)=0,\!5403;\\[2pt] & p_{0}^{0}= p(0;0)=1; & & p_{1}^{0}= \nu(0,\!5;0)=2; & & p_{2}^{0}= p(1;0)=3. \end{aligned}$

2. Выведем разностные уравнения, определяющие положение точек на следующем временном слое и значения искомых функций в них. Система (8.74) заменяется системой разностных уравнений. Численное решение задачи строится на характеристической сетке, которая формируется одновременно с процессом нахождения решения. При этом обычно рассматривается трехточечный угловой шаблон с узлами "1","2","3", составленный отрезками характеристик (из семейств $C^{+}$ и $C^{-}$ ), пересекающихся в рассчитываемой точке "3" (рис. 8.19,б). Назовем этот шаблон базовым. Предполагается, что положение узлов "1" , "2" и значения решения в узлах известны из начальных условий или рассчитаны на предыдущем шаге. Положение узла "3" и значение решения в нем подлежат определению.

Аппроксимируя производную по формулам

${\left(\frac{dx}{dt}\right)}_{1}\cong \frac{x_3-x_1}{t_3-t_1}$ (вдоль характеристики из семейства $C^{+}$ с $\lambda_1=1$ ),

${\left(\frac{dx}{dt}\right)}_{2}\cong \frac{x_3-x_2}{t_3-t_2}$ (вдоль характеристики из семейства $C^{-}$ с $\lambda_2=-1$ ),

получаем

$x_3-x_1=t_3-t_1,\qquad x_3-x_2=t_2-t_3.$

(8.75)

Условие $p-\nu=\text{const}$ вдоль характеристик из семейств $C^{+}$ записывается для лежащих на ней узлов "1" и "3":

$p_3-\nu_3=p_1-\nu_1,$

(8.76)

а условие $p-\nu=\text{const}$ вдоль характеристики из семейства $C^{-}$ записывается для лежащих на ней узлов "2" и "3":

$p_3+\nu_3= p_2+\nu_2.$

(8.77)

Система (8.75)–(8.77) является замкнутой относительно искомых значений $x_3,\,t_3,\,p_3,\,\nu_3$ при известных $x_1, t_1, x_2, t_2, p_1, \nu_1, p_2, \nu_2$ . Формулы (8.75)–(8.77) являются базовыми для расчета всех точек, принадлежащих области $D$ . После совмещения узлов "1" и "2" базового шаблона с выбранными очередными узлами сетки в (8.75)–(8.77) подставляются конкретные величины в соответствии с принятыми обозначениями.

3. Производится расчет первого временного слоя. В решаемой задаче на нулевом временном слое угловой шаблон может занимать два положения (рис. 8.20, а,б). Для точки с номерами $(0,1)$ первого рассчитываемого слоя узлы "1" и "2" углового шаблона совмещаются с узлами $(0,0)$ и $(1,0)$ (рис. 8,20,а). Поэтому система (8.75)–(8.77) имеет вид

$\begin{aligned} & x_{0,1}-x_{0,0}= t_{0,1}-t_{0,0}, & & x_{0,1}-x_{1,0}= t_{1,0}-t_{0,1},\\[2pt] & p_0^1-\nu_0^1= p_0^0-\nu_0^0, & & p_0^1+ \nu_0^1= p_1^0+ \nu_1^0. \end{aligned}$

Подставляя $t_{0,0}=0,~ t_{1,0}=0,~ x_{0,0}=0,~ x_{1,0}=0,\!5$ и известные значения $p_0^0, p_1^0, \nu_0^0, \nu_1^0$ на нулевом слое (см. п.1), получаем

$\begin{gathered}x_{0,1}-0=t_{0,1}-0;\qquad x_{0,1}-0,\!5=0-t_{0,1};\\[2pt] p_0^1-\nu_0^1=1-1=0;\qquad p_0^1+\nu_0^1=2+ 0,\!8776= 2,\!8776. \end{gathered}$

Решение этой системы: $x_{0,1}=0,\!25;~ t_{0,1}=0,\!25;~ p_0^1=1,\!4388;~ \nu_0^1=1,\!4388$ .

Для второй точки первого слоя (рис. 8.20,б) с номерами $(l;l)$ узлы "1" и "2" углового шаблона совмещаются с узлами $(l;0)$ и $(2;0)$ . Тогда система (8.75)–(8.77) имеет вид:

$\begin{gathered}x_{1,1}-x_{1,0}=t_{1,1}-t_{1,0};\qquad x_{1,1}-x_{2,0}= t_{2,0}-t_{1,1};\\[2pt] p_1^1-\nu_1^1= p_1^0-\nu_1^0;\qquad p_1^1+ \nu_1^1= p_2^0+ \nu_2^0. \end{gathered}$

Подставляя $t_{1,0}=0,~ t_{2,0}=0,~ x_{1,0}=0,\!5,~ x_{2,0}=1$ и известные значения $p_1^0,\nu_1^0,o_2^0,\nu_2^0$ на нулевом слое (см. п.1), получаем

$\begin{gathered}x_{1,1}-0,\!5=t_{11}-0;\qquad x_{1,1}-1=0-t_{1,1};\\[2pt] p_1^1-\nu_1^1=2-0,\!8776=1,\!1224;\qquad p_1^1+\nu_1^1=3+0,\!5403= 3,\!5403. \end{gathered}$

Отсюда $x_{1,1}= 0,\!75;~ t_{1,1}= 0,\!25; p_1^1= 2,\!33135;~ \nu_1^1= 1,\!20895$ .

4. Для расчета точки $(0,2)$ следующего временного слоя узлы "1" и "2" углового шаблона совмещаются с узлами $(0,1)$ и $(1,1)$ (рис.8.20,в). Тогда система (8.75)–(8.77) принимает вид:

$\begin{gathered}x_{0,2}-x_{0,1}=t_{0,2}-t_{0,1};\qquad x_{0,2}-x_{1,1}=t_{1,1}-t_{0,2};\\[2pt] p_0^2-\nu_0^2= p_0^1-\nu_0^1;\qquad p_0^2+\nu_0^2= p_1^1+ \nu_1^1. \end{gathered}$

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

$\begin{gathered}x_{0,2}-0,\!25=t_{0,2}-0,\!25;\qquad x_{0,2}-0,\!75= 0,\!25-t_{0,2};\\[2pt] p_0^2-\nu_0^2= 1,\!4388-1,\!4388;\qquad p_0^2+ \nu_0^2=2,\!33135+ 1,\!20895. \end{gathered}$

Отсюда $x_{0,2}= 0,\!5;~ t_{0,2}= 0,\!5;~ p_0^2=\nu_0^2= 1,\!77015$ .

Очевидно (см. рис. 8.19,а), количество узлов на каждом последующем слое уменьшается на единицу. Расчет продолжается до замыкания области, так как на последнем слое останется один узел. В данном примере этим узлом является точка $(0,2)$ .

Замечание. Рассмотренная в примере 8.10 задача является всего лишь модельной, выбранной для пояснения основных принципов конструирования метода характеристик и его численной реализации. В вычислительной практике, как правило, ставятся нелинейные задачи. Проиллюстрируем некоторые особенности их решения на следующем примере.

Пример 8.11. Получить соотношения метода характеристик для решения нелинейной задачи Коши

$\begin{gathered}\begin{aligned}&a_{11}\frac{\partial u_1(x,t)}{\partial t}+ a_{12}\frac{\partial u_2(x,t)}{\partial t}+ b_{11}\frac{\partial u_1(x,t)}{\partial x}+ b_{12}\frac{\partial u_2(x,t)}{\partial x}=f_1,\\ &a_{21}\frac{\partial u_1(x,t)}{\partial t}+ a_{22}\frac{\partial u_2(x,t)}{\partial t}+ b_{21}\frac{\partial u_1(x,t)}{\partial x}+ b_{22}\frac{\partial u_2(x,t)}{\partial x}=f_2,\end{aligned}\quad (x,t)\in D\\[4pt] u_1(x,0)= \varphi_1(x),\quad u_2(x,0)= \varphi_2(x),\quad x\in \mathbb{R},\end{gathered}$

где $a_{ij}=a_{ij}(x,t,u),~ b_{ij}=b_{ij}(x,t,u),~ f_{i}=f_{i}(x,t,u),~ i=1,2;~ j=1,2;~ \varphi_{1}(x),~ \varphi_2(x)$ — заданные функции; $D=\mathbb{R}\times (0,T),~ T>0$ — заданный момент времени.

[reshenie=Решение, формулы (8.78)-(8.80)]

Поставленная задача является частным случаем задачи (8.65) при $m=2$ . Применим метод характеристик, включающий два этапа.

Первый этап. Нахождение характеристик и характеристических соотношений. Уравнение (8.67) имеет вид

$|B-A \lambda|= \begin{vmatrix}b_{11}-\lambda a_{11}& b_{12}-\lambda a_{12}\\ b_{21}-\lambda a_{21}& b_{22}-\lambda a_{22}\end{vmatrix}=0$ или $a\,\lambda^2+ b\,\lambda+c=0,$

где $a=a_{11}a_{12}-a_{12}a_{21},~ b=a_{21}b_{12}+ a_{12}b_{21}-a_{11}b_{22}-a_{22}b_{11},~ c=b_{11}b_{22}-b_{12}b_{21}$ . Решая полученное уравнение, находим выражения для двух характеристических направлений в каждой точке расчетной области:

$\begin{aligned}&{\left(\frac{dx}{dt}\right)}_{1}= \lambda_{1}(x,t,u_1,u_2)= \frac{-b+ \sqrt{b^2-ac}}{2a}\quad (C^{+});\\[2pt] &{\left(\frac{dx}{dt}\right)}_{2}= \lambda_{2}(x,t,u_1,u_2)= \frac{-b-\sqrt{b^2-ac}}{2a}\quad (C^{-}). \end{aligned}$

(8.78)

Сформируем расширенную матрицу

$\bigl(B\,dt-A\,dx,\, f\,dt-A\,du\bigr)= \begin{pmatrix}b_{11}dt-a_{11}dx& b_{12}dt-a_{12}dx& f_1dt-a_{11}du_1-a_{12}du_2\\ b_{21}dt-a_{21}dx& b_{22}dt-a_{22}dx& f_2dt-a_{21}du_1-a_{22}du_2 \end{pmatrix}$

Характеристические соотношения получаются в результате приравнивания нулю определителя второго порядка $(m=2)$ , составленного из второго и третьего столбцов расширенной матрицы:

$\begin{vmatrix}b_{12}dt-a_{12}dx& f_1dt-a_{11}du_1-a_{12}du_2\\ b_{22}dt-a_{22}dx& f_2dt-a_{21}du_1-a_{22}du_2 \end{vmatrix}=0.$

Раскрыв определитель и поделив на $dt\ne0$ , получим

$\left(e-a\,\frac{dx}{dt}\right)\!du_1+ g\,du_2+ h\,dt+ q\,dx=0,$

где $e=a_{11}b_{22}-a_{21}b_{12},~ g=a_{12}b_{22}-a_{22}b_{12},~ h=b_{12}f_2-b_{22}f_1,~ q=a_{22}f_1-a_{12}f_2$ . Последние два слагаемых перепишем в эквивалентной форме $h\,dt+ q\,dx= \left(h+ q\,\frac{dx}{dt}\right)\!dt$ .

Поскольку для семейства характеристик $C^{+}$ выполняется равенство ${\left(\frac{dx}{dt}\right)}_{1}= \lambda_{1}$ , а для семейства $C^{-}$ — равенство ${\left(\frac{dx}{dt}\right)}_{2}= \lambda_{2}$ , то характеристические \соотношения имеют вид

$\begin{aligned}& (e-a \lambda_1)\,du_1+ g\,du_2+ (h+q \lambda_1)\,dt=0\quad (C^{+});\\[2pt] & (e-a \lambda_2)\,du_1+ g\,du_2+ (h+q \lambda_2)\,dt=0\quad (C^{-}). \end{aligned}$

(8.79)

В результате получена характеристическая система (8.78), (8.79) из четырех обыкновенных дифференциальных уравнений. Они позволяют найти значения $t,x,u_1,u_2$ вдоль характеристик.

Второй этап. Конструирование численной схемы. Как и в предыдущем примере, будем использовать трехточечный угловой шаблон (рис. 8.19, б). Предполагается, что в точках "1" и "2" известны значения $t,x,u_1,u_2$ из начальных данных или из расчета предыдущего слоя, т.е. $t_1,x_1,u_{11},u_{21}$ для точки "1" и $t_2,x_2,u_{12},u_{22}$ для точки "2". Здесь при обозначении значений функций $u_1,u_2$ , а далее и при обозначении других функций в точках шаблона первый нижний индекс соответствует номеру самой функции, а второй индекс — номеру точки в шаблоне. Требуется найти значения $t_3,x_3, u_{13}, u_{23}$ определяющие положение точки "3" и значения искомых функций в ней. Предполагается, что точки "1" и "2" близки друг к другу, поэтому заменим дифференциалы конечными приращениями, а производную $\frac{dx}{dt}$ — конечноразностным отношением:

$\begin{gathered}dt\cong \Delta t=t_3-t_1,\quad dx\cong \Delta x=x_3-x_1,\quad du_1\cong \Delta u_1=u_{13}-u_{11},\quad du_2\cong \Delta u_2=u_{23}-u_{21},\quad {\left(\frac{dx}{dt} \right)}_{1}\cong \frac{x_3-x_1}{t_3-t_1}\quad (C^{+});\\[4pt] dt\cong \Delta t=t_3-t_2,\quad dx\cong \Delta x=x_3-x_2,\quad du_1\cong \Delta u_1=u_{13}-u_{12},\quad du_2\cong \Delta u_2=u_{23}-u_{22},\quad {\left(\frac{dx}{dt}\right)}_{2}\cong \frac{x_3-x_2}{t_3-t_2}\quad (C^{+}). \end{gathered}$

В итоге уравнения характеристик являются уравнениями прямых, соединяющих точки "1" и "2" с точкой "3", а характеристическая система (8.78), (8.79) преобразуется к виду

$\begin{gathered}\begin{aligned}& x_3-x_1= (\lambda_1)_1(t_3-t_1),\\ & x_3-x_2= (\lambda_2)_2(t_3-t_2),\end{aligned}\\[4pt] \begin{aligned}& (e-a \lambda_1)_1(u_{13}-u_{11})+ (g)_1(u_{23}-u_{21})+ (h+ q \lambda_1)_1 (t_3-t_1)=0,\\ & (e-a \lambda_2)_2(u_{13}-u_{12})+ (g)_2(u_{23}-u_{22})+ (h+ q \lambda_2)_2 (t_3-t_2)=0. \end{aligned}\end{gathered}$

(8.80)

Здесь нижние индексы 1 и 2 указывают, что значения всех функций, входящих в соответствующие выражения, подсчитываются в точке "1" или "2". Решая систему (8.80), можно найти значения $t_3,x_3,u_{13},u_{23}$ " определяющие положение точки "3" и значения искомых функций в ней. Сначала обычно решаются первые два уравнения и определяются значения $t_3,x_3$ , а затем — последние два уравнения, откуда находятся $u_{13},u_{23}$ .

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

1. При $t=0$ на оси абсцисс задаются точки с некоторым шагом. В этих точках (на рис. 8.21 выделен фрагмент, содержащий точки 1-6) значения функций $u_{1},u_{2}$ определяются начальными данными, т.е. функциями $\varphi_1(x), \varphi_2(x)$ .

Схема численного решения уравнения в частных производных

2. В соответствии с изложенной методикой последовательно рассчитываются координаты точек 7-11 и приближенные значения решения в них. При этом точки "1" и "2" углового шаблона последовательно совмещаются с точками 1 и 2, затем с точками 2 и 3, 3 и 4, 4 и 5, 5 и 6. Узлами выстраиваемой характеристической сетки служат точки пересечения характеристик (они образуют следующий временной слой).

3. Расчет слоев продолжается до достижения последней точки (точки 21) треугольной области.

Отметим, что дифференциальным уравнениям (8.78) соответствуют криволинейные характеристики, а точка "3" должна быть точкой пересечения этих криволинейных характеристик. При переходе от дифференциальных уравнений (8.78),(8.79) к разностным дифференциалы были заменены конечными разностями, а производные конечно-разностными отношениями. Вследствие этого криволинейные характеристики были аппроксимированы отрезками прямых. Решение системы (8.80) объявляется первым приближением, т.е. $t_3^{(1)}= t_3,~ x_3^{(1)}= x_3,$ $u_{13}^{(1)}= u_{13},~ u_{23}^{(1)}= u_{23},$ и может возникнуть необходимость в уточнении координат точки "3" и значений искомых функций в этой точке. Ограничимся изложением одного из возможных способов уточнения. После нахождения первого приближения следует найти средние арифметические значения:

$\begin{aligned}& \overline{\lambda}_{1}= \frac{(\lambda_1)_1+ (\lambda_1)_3}{2}\,,\quad \overline{\lambda}_{2}= \frac{(\lambda_2)_2+ (\lambda_2)_3}{2}\,,\quad K_1= \frac{(e-a \lambda_1)_1+ (e-a \lambda_1)_3}{2}\,,\quad K_2= \frac{(e-a \lambda_2)_2+ (e-a \lambda_2)_3}{2}\,,\\[4pt] & \overline{L}_{1}= \frac{(g)_1+ (g)_3}{2}\,,\quad \overline{L}_{2}= \frac{(g)_2+ (g)_3}{2}\,,\quad M_1= \frac{(h+ q \lambda_1)_1+ (h+ q \lambda_1)_3}{2}\,,\quad M_2= \frac{(h+ q \lambda_2)_2+ (h+ q \lambda_2)_3}{2}\,,\end{aligned}$

где нижний индекс 3 указывает, что соответствующее выражение подсчитывается в найденной точке "3". Эти соотношения получаются, если при интегрировании дифференциальных уравнений применять правило трапеций. Затем решается система

$\begin{gathered}x_{3}^{(2)}-x_{1}= \overline{\lambda}_{1}(t_{3}^{(2)}-t_{1}),\\ x_{3}^{(2)}-x_{2}= \overline{\lambda}_{2}(t_{3}^{(2)}-t_{2}),\\[4pt] K_1\,\bigl(u_{13}^{(2)}-u_{11}\bigr)+ L_1\,\bigl(u_{23}^{(2)}-u_{21}\bigr)+ M_1\,\bigl(t_{3}^{(2)}-t_{1}\bigr)=0,\\ K_2\,\bigl(u_{13}^{(2)}-u_{12}\bigr)+ L_2\,\bigl(u_{23}^{(2)}-u_{22}\bigr)+ M_2\,\bigl(t_{3}^{(2)}-t_{2}\bigr)=0 \end{gathered}$

относительно нового (второго) приближения $t_2^{(2)},\,x_3^{(2)},\, u_{13}^{(2)},\, u_{23}^{(2)}$ . Описанный процесс продолжают до тех пор, пока решения, полученные при двух последовательных приближениях, не будут совпадать с заданной точностью. Как правило, если расстояние между точками "1" и "2" невелико, достаточно сделать два уточнения.

Замечания

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

2. Описанный метод расчета, когда точки следующего слоя заранее неизвестны, а вычисления ведутся вдоль характеристик семейств $C^{+}$ и $C^{-}$ , выделяемых в процессе расчета, называется классическим методом характеристик. Существуют модификации этого метода, когда точки следующего слоя задаются заранее, а для нахождения значений искомых функций в них характеристики выпускаются в направлении к известному слою. Такие методы называются сеточно-характеристическими.

Метод Годунова

В механике сплошных сред при решении задач, описывающихся квазилинейными дифференциальными уравнениями в частных производных, в расчетных областях могут возникать поверхности сильного разрыва, к которым относятся ударные волны, центрированные волны разрежения, поверхности контактного разрыва. На этих поверхностях терпят разрыв как искомые функции, характеризующие исследуемые процессы, так и их производные. При решении таких задач используются два класса численных схем, в одном из которых с использованием законов сохранения явно выделяются поверхности разрыва, а в другом эти поверхности явно не выделяются, а рассчитываются "насквозь" вместе с областями непрерывности. К последнему классу методов относится метод С. К. Годунова. Идею его конструирования рассмотрим на примере решения смешанной краевой задачи для уравнений акустики.

Пример 8.12. Получить соотношения метода Годунова для решения задачи

$\frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0$

с начальными условиями $u(x,0)=\psi_1(x),~ p(x,0)=\psi_2(x),~ 0 \leqslant x \leqslant 1$ ,
и краевыми условиями $p(0,t)=\varphi_1(t),~ p(1,t)=\varphi_2(t),~ 0<t \leqslant T$ .

Решение, формулы (8.81)-(8.90)

Система двух дифференциальных уравнений описывает распространение плоских звуковых волн (малых возмущений) в покоящейся среде. Здесь $u(x,t),\, p(x,t)$ — малые отклонения скорости и давления от их значений в невозмущенной среде, вызванные распространением звуковых волн; $\rho_0$ — плотность невозмущенной среды; $a_0$ — скорость звука.

Рассмотрим основные этапы конструирования метода Годунова.

1. Введем разбиение отрезка $[0;1]$ с шагом $h=1\!\!\not{\phantom{|}}\,n$ , где $n$ — число промежутков. В результате получим узлы $x_{i},~ i=0,1,\ldots,n$ , определяемые соотношением

$x_{i+1}=x_{i}+h,\quad i=0,1,\ldots,n-1;\quad x_0=0,\quad x_n=1.$

2. Заменим функции $\psi_1(x),\, \psi_2(x)$ на начальном слое $t=0$ кусочно-постоянными (пример описанной замены для некоторой произвольной функции $\psi_1(x)$ изображен на рис. 8.22) путем их интегрального осреднения по формулам

$u_{1-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} \psi_1(x)\,dx\,,\quad p_{1-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} \psi_2(x)\,dx\,,\quad i=1,2,\ldots,n.$

(8.81)

Здесь и далее все заданные или ранее рассчитанные функции снабжаются нижними полуцелыми индексами. Из рис. 8.22 видно, что полученная кусочно-постоянная функция имеет в узлах разрывы первого рода.

3. Выберем шаг $\tau$ по времени $t$ и образуем следующий расчетный слой $(t_1=\tau)$ с теми же узлами $x_{i}$ (рис. 8.23). Таким образом, образуется совокупность прямоугольных ячеек. Предполагается, что искомые функции на горизонтальных и вертикальных границах образованных ячеек также являются постоянными.

На рассчитываемом слое при $t=\tau$ приближенное решение задачи также представляется в виде кусочно-постоянной функции, имеющей разрывы в узлах сетки, совпадающих с узлами на предыдущем слое. Это означает, что во времени на слоях $t=\text{const}$ сохраняется структура решения, изображенная на рис. 8.22. Значения искомых функций на рассчитываемом слое (верхних границах ячеек) обозначаются верхними полуцелыми индексами: $u^{1-1\!\not{\phantom{|}}\,\,2},\, p^{1-1\!\not{\phantom{|}}\,\,2},$ $i=1,2,\ldots,n$ . На вертикальных границах ячеек значения функций обозначаются прописными буквами с нижними целыми индексами: $U_{i},\,P_{i},~ i=0,\ldots,n$ .

При $i=0,~ i=n$ значения $P_0$ и $P_n$ определяются краевыми условиями по формулам, аналогичным (8.81):

$P_0= \frac{1}{\tau} \int\limits_{0}^{\tau} \varphi_1(t)\,dt\,,\quad P_n= \frac{1}{\tau} \int\limits_{0}^{\tau} \varphi_2(t)\,dt\,.$

(8.82)

Значения $U_1,\ldots,U_{n-1},\, P_1,\ldots,P_{n-1}$ находятся в результате решения вспомогательной задачи о распаде произвольного разрыва, описанного в п.4, а значения $U_0,\,U_n$ определяются в п.5.

4. Выберем две соседние ячейки с узлами $x_{i-1},\,x_{i}$ и $x_{i},\,x_{i+1}$ (эти ячейки выделены на рис. 8.23 жирными линиями). С целью нахождения значений $U_{i},\,P_{i}$ на внутренней границе ячеек при $x=x_{i}$ рассмотрим решение вспомогательной задачи, называемой задачей о распаде произвольного разрыва.

Дана система уравнений

$\frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0$

(8.83)

с начальными условиями вида

$u(x,0)= \left\{\begin{aligned}&u_{i-1\!\not{\phantom{|}}\,\,2},\quad x_{i-1} \leqslant x<x_{i},\\ &u_{i+1\!\not{\phantom{|}}\,\,2},\quad x_{i} \leqslant x<x_{i+1}, \end{aligned}\right.\qquad p(x,0)= \left\{\begin{aligned}&p_{i-1\!\not{\phantom{|}}\,\,2},\quad x_{i-1} \leqslant x<x_{i},\\ &p_{i+1\!\not{\phantom{|}}\,\,2},\quad x_{i} \leqslant x<x_{i+1}. \end{aligned}\right.$

Требуется найти решение $u(x,t),\,p(x,t)$ при $t>0$ и $x=x_{i}$ , то есть $u(x_{i},t)=U_{i},~ p(x_{i},t)=P_{i}$ .

Заметим, что в данной задаче начальные данные имеют разрыв первого рода в точке $x_{i}$ . При $t>0$ они трансформируются, т.е. происходит распад разрыва, что обусловливает название рассматриваемой задачи. Для ее решения применим метод характеристик. Перепишем систему (8.83) в матричной форме

$A\cdot \frac{\partial u(x,t)}{\partial t}+ B\cdot \frac{\partial u(x,t)}{\partial x}=0$ , где $A=\begin{pmatrix}1&0\\ 0&1 \end{pmatrix}\!,~~ B=\begin{pmatrix}0& 1\!\!\not{\phantom{|}}\,\rho_0\\ \rho_0a_0^2&0 \end{pmatrix}\!,~~ u=\begin{pmatrix}u\\p \end{pmatrix}$ .

Для нахождения характеристических направлений используем условие

$|B-A\cdot \lambda|= \begin{vmatrix}-\lambda& 1\!\!\not{\phantom{|}}\,\rho_0\\ \rho_0a_0^2&-\lambda\end{vmatrix}=0$ . Отсюда ${\left(\frac{dx}{dt}\right)\!}^{+}= \lambda^{+}= a_0,~~ {\left(\frac{dx}{dt}\right)\!}^{-}= \lambda^{-}=-a_0$ .

Интегрируя последние уравнения, имеем

$x-a_0t=\text{const}~(C^{+}),\qquad x+a_0t=\text{const}~(C^{-}).$

Если при $t=0$ выполняется условие $x=x_{i}$ , то из последних соотношений выделяются две прямые, исходящие из точки $x_{i}$ на начальном слое:

$x-a_0t=x_{i}~(C^{+}),\qquad x+a_0t=x_{i}~(C^{-}).$

Для получения характеристических соотношений запишем расширенную матрицу (8.69) (учитывая, что $f\equiv 0$ ):

$\left(B\,dt-A\,dx,-A\! \begin{pmatrix}du\\dp\end{pmatrix}\right)= \begin{pmatrix}-dx&\dfrac{1}{\rho_0}\,dt&-dt\\ \rho_0a_0^2dt&-dx&-dp \end{pmatrix}\!.$

Приравнивая определитель матрицы, составленной из двух последних столбцов, нулю, получаем

$-\dfrac{1}{\rho_0}\,dp\,dt-dx\,du=0.$ Полагая $dt\ne0$ , разделим уравнение на $dt\colon\,-\dfrac{1}{\rho_0}\,dp-\frac{dx}{dt}\,du=0.$

Подставляя вместо $\frac{dx}{dt}$ сначала $\lambda^{+}=a_0$ , а затем $\lambda^{-}=-a_0$ , получаем

$du+ \dfrac{1}{\rho_0\,a_0}\,dp=0\quad (C^{+}),\qquad du-\dfrac{1}{\rho_0\,a_0}\,dp= 0\quad (C^{-}).$

Так как коэффициенты обоих уравнений постоянны, то их можно переписать в эквивалентной форме

$d\! \left(u+ \dfrac{1}{\rho_0\,a_0}\,p\right)=0\quad (C^{+}),\qquad d\! \left(u-\dfrac{1}{\rho_0\,a_0}\,p\right)=0\quad (C^{-}).$

Интегрируя и обозначая постоянные интегрирования $I^{+},\,I^{-}$ , имеем

$u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}\quad (C^{+}),\qquad u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}\quad (C^{-}).$

(8.84)

Полученные соотношения называются инвариантами Римана. Они показывают, что величина $u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}$ остается постоянной вдоль семейства $C^{+}$ , а величина $u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}$ остается постоянной вдоль семейства $C^{-}$ .

Теперь, используя направления характеристик и сохраняющиеся вдоль них инварианты, построим численную схему решения задачи о распаде разрыва.

Семества прямых, исходящие из узлов сетки

Рассмотрим две характеристики $x-a_0t=x_{i-1}$ и $x-a_0t=x_{i}$ из семейства $C^{+}$ , исходящие из узлов $x_{i-1},\,x_{i}$ , а также две характеристики $x+a_0t=x_{i}$ и $x+a_0t=x_{i+1}$ из семейства $C^{-}$ , исходящие из узлов $x_{i},\,x_{i+1}$ (рис. 8.24). В результате пересечения указанных характеристик образуются области I, II, III , первые две из которых представлены треугольниками, а третья — четырехугольником. Вертикальная граница $(x_{i}= \text{const})$ , разделяющая ячейки, на которой ищется решение задачи, располагается в области III. Покажем, что на этой границе искомое решение постоянно. Для этого выберем на ней произвольную точку (назовем ее точка 3) , а решение в ней обозначим $U_{i},\, P_{i}$ . В эту точку приходит характеристика семейства $C^{+}$ , исходящая из точки 1 промежутка $(x_{i-1}, x_{i})$ . Точке 1 соответствуют заданные значения $u_{i-1\!\not{\phantom{|}}\,\,2},\, p_{i-1\!\not{\phantom{|}}\,\,2}$ . Вдоль этой характеристики остается постоянной величина $u+ \dfrac{1}{\rho_0\,a_0}\,p=I^{+}$ , то есть

$U_{i}+ \frac{1}{\rho_0\,a_0}\cdot P_{i}= u_{i-1\!\not{\phantom{|}}\,\,2}+ \frac{1}{\rho_0\,a_0}\cdot p_{i-1\!\not{\phantom{|}}\,\,2}.$

Кроме того, в точку 3 приходит характеристика семейства $C^{-}$ , исходящая из точки 2 промежутка $(x_{i},x_{i+1})$ . Точке 2 соответствуют заданные значения и $u_{i+ 1\!\not{\phantom{|}}\,\,2},\, p_{i+ 1\!\not{\phantom{|}}\,\,2}$ . Вдоль этой характеристики остается постоянной величина $u-\dfrac{1}{\rho_0\,a_0}\,p=I^{-}$ , то есть

$U_{i}-\frac{1}{\rho_0\,a_0}\cdot P_{i}= u_{i+ 1\!\not{\phantom{|}}\,\,2}-\frac{1}{\rho_0\,a_0}\cdot p_{i+ 1\!\not{\phantom{|}}\,\,2}.$

Разрешая систему из двух последних уравнений относительно $U_{i}$ и $P_{i}$ , получаем

$\begin{gathered}U_{i}= \frac{u_{i-1\!\not{\phantom{|}}\,\,2}+ u_{i+ 1\!\not{\phantom{|}}\,\,2}}{2}-\frac{p_{i+ 1\!\not{\phantom{|}}\,\,2}-p_{i+ 1\!\not{\phantom{|}}\,\,2}}{2\rho_0\,a_0}\,\\ P_{i}= \frac{p_{i-1\!\not{\phantom{|}}\,\,2}+ p_{i+ 1\!\not{\phantom{|}}\,\,2}}{2}+ \rho_0\,a_0\cdot \frac{u_{i+ 1\!\not{\phantom{|}}\,\,2}-u_{i-1\!\not{\phantom{|}}\,\,2}}{2}\,,\end{gathered}\quad i=1,2,\ldots,n-1.$

(8.85)

Поскольку точка 3 была выбрана произвольно на вертикальной границе ячейки, то полученные выражения определяют решение в любой точке этой границы (очевидно, решение постоянно вдоль границы), за исключением самой точки $x_{i}$ и точки пересечения характеристик, исходящих из точек $x_{i-1},\, x_{i+1}$ . Таким образом, формулы (8.85) определяют "большие" величины $U_{i},\, P_{i}$ на всех внутренних вертикальных границах ячеек.

Результаты решения данной вспомогательной задачи определяют шаг $\tau$ по времени. Действительно, влияние каждого из разрывов, относящихся к соседним узлам $x_{i-1},\,x_{i}$ на слое $t=0$ , ограничивается характеристиками, исходящими из этих узлов. Физически это означает, что разрывы, возникшие при $t=0$ , распространяются со скоростью звука $a_0$ . Характеристики, исходящие из соседних узлов навстречу друг другу, пересекутся в точке с координатами $t=\frac{h}{2a_0},~ x=\frac{x_{i-1}+ x_{i}}{2}$ . От этой встречи разрывов возникнут новые волны, которые через следующий отрезок времени длительностью $\frac{h}{2a_0}$ достигнут вертикальных границ $(x_{i-1}=\text{const},~ x_{i}=\text{const})$ соседних ячеек. Именно до этого момента, т.е. в течение времени $\frac{h}{a_0}= \frac{h}{2a_0}+ \frac{h}{2a_0}$ , значения $U_{i},\, P_{i}$ будут сохраняться неизменными и определяться по формулам (8.85). Таким образом, величина шага по времени должна выбираться из условия $\tau< \frac{h}{a_0}$ .

При $i=0$ для нахождения величины $U_0$ при известном $P_0$ используется инвариант Римана $u-\frac{p}{\rho_0 a_0}=I^{-}$ вдоль семейства $C^{-}$ (рис. 8.25). Поскольку на промежутке $(0,x_1)$ заданы значения $u_{1\!\not{\phantom{|}}\,\,2},\,p_{1\!\not{\phantom{|}}\,\,2}$ , то $U_0-\frac{P_0}{\rho_0a_0}= u_{1\!\not{\phantom{|}}\,\,2}-\frac{p_{1\!\not{\phantom{|}}\,\,2}}{\rho_0 a_0}$ . Отсюда получаем

$U_0= \frac{1}{\rho_0a_0}\cdot P_0+ u_{1\!\not{\phantom{|}}\,\,2}-\frac{1}{\rho_0 a_0}\cdot p_{1\!\not{\phantom{|}}\,\,2}.$

(8.86)

Семейства прямых, образующих прямоугольные области

При $i=n$ для нахождения величины $U_n$ при известном $P_n$ используется инвариант Римана $u+\frac{p}{\rho_0a_0}=I^{+}$ вдоль семейства $C^{+}$ (рис. 8.25). Поскольку на последнем промежутке $(x_{n-1},1)$ заданы значения $u_{n-1\!\not{\phantom{|}}\,\,2},\, p_{n-1\!\not{\phantom{|}}\,\,2}$ , то $U_n+\frac{P_n}{\rho_0a_0}= u_{n-1\!\not{\phantom{|}}\,\,2}+ \frac{P_{n-1\!\not{\phantom{|}}\,\,2}}{\rho_0a_0}$ . Отсюда получаем

$U_n=-\frac{1}{\rho_0a_0}\cdot P_n+ u_{n-1\!\not{\phantom{|}}\,\,2}+ \frac{1}{\rho_0a_0}\cdot p_{n-1\!\not{\phantom{|}}\,\,2}.$

(8.87)

На рис. 8.25 показаны две крайние ячейки с тремя характеристиками, одна из которых (верхняя) ограничивает область влияния начальных данных, заданных на промежутках $(x_0,x_1)$ и $(x_{n-1},x_{n})$ соответственно.

6. Перейдем к задаче расчета искомых величин на следующем временном слое. Для этого рассмотрим одну ячейку (рис. 8.26), для которой известны значения приближенного решения $u_{i-1\!\not{\phantom{|}}\,\,2},\, p_{i-1\!\not{\phantom{|}}\,\,2}$ на нижней границе и $U_{i-1},P_{i-1},U_{i},P_{i},$ $i=1,2,\ldots,n$ , на левой и правой вертикальных границах (см. п.4). Требуется определить значения $p^{i-1\!\not{\phantom{|}}\,\,2},\, u^{i-1\!\not{\phantom{|}}\,\,2}$ на верхней границе ячейки.

Для решения поставленной задачи перепишем уравнения (8.83) в так называемом дивергентном виде

$\frac{\partial u(x,t)}{\partial t}+ \frac{\partial}{\partial x}\! \left(\frac{p(x,t)}{\rho_0} \right)=0,\qquad \frac{\partial p(x,t)}{\partial t}+ \frac{\partial}{\partial x}\bigl(\rho_0a_0^2u(x,t) \bigr)=0.$

(8.88)

Запишем двойной интеграл от левой и правой частей уравнений по поверхности ячейки. Для его нахождения применим формулу Грина перехода от двойного интеграла к криволинейному интегралу второго рода:

$\iint\limits_{S}\! \left(\frac{\partial \Phi_1}{\partial t}+ \frac{\partial \Phi_2}{\partial x}\right)\!dx\,dt= \oint\limits_{\gamma}\bigl(\Phi_1\,dx-\Phi_2\,dt\bigr),$

где $S$ — поверхность ячейки, а $\gamma$ — ее граница, обходимая против часовой стрелки. Таким образом, вместо системы (8.88) получается система интегральных уравнений

$\oint\limits_{\gamma}u\,dx-\frac{p}{\rho_0}\,dt=0,\qquad \oint\limits_{\gamma}p\,dx-\rho_0a_0^2u\,dt=0.$

Проводя интегрирование, имеем

$\begin{gathered}u_{i-1\!\not{\phantom{|}}\,\,2}h-\frac{P_{i}}{\rho_0}\cdot \tau-u^{i-1\!\not{\phantom{|}}\,\,2}h+ \frac{P_{i-1}}{\rho_0}\cdot \tau=0,\\[4pt] p_{i-1\!\not{\phantom{|}}\,\,2}h-\rho_0a_0^2U_{i}\tau-p^{i-1\!\not{\phantom{|}}\,\,2}h+ \rho_0a_0^2U_{i-1}\tau=0. \end{gathered}$

Из этих уравнений получаем

$\begin{gathered}u^{i-1\!\not{\phantom{|}}\,\,2}= u_{i-1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}\cdot (P_{i}-P_{i-1}),\\[4pt] p^{i-1\!\not{\phantom{|}}\,\,2}= p_{i-1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\cdot \rho_0a_0^2(U_{i}-U_{i-1}). \end{gathered}$

Таким образом, получены формулы для расчета первого временного слоя $(t=\tau)$ . После нахождения значений $u^{\frac{1}{2}},p^{\frac{1}{2}}, u^{\frac{3}{2}}, p^{\frac{3}{2}},\ldots, u^{n-\frac{1}{2}}, p^{n-\frac{1}{2}}$ из условия $\tau<\frac{h}{a_0}$ выбирается очередной шаг и формируется новый временной слой $(t_2=t_1+ \tau)$ . Найденные значения используются в качестве опорных (они переобозначаются путем присвоения им соответствующих нижних индексов) для нахождения приближенного решения на следующем слое по формулам (8.85)–(8.87), (8.89). Заметим, что для вычисления $P_0,\,P_n$ вместо (8.82) используются формулы интегрального осреднения, аналогичные (8.81). Для расчета, например, второго временного слоя они имеют вид

$P_0= \frac{1}{\tau} \int\limits_{t_1}^{t_2} \varphi_1(t)\,dt\,,\qquad P_n= \frac{1}{\tau} \int\limits_{t_1}^{t_2} \varphi_2(t)\,dt\,.$

(8.90)

Изложенная схема метода Годунова имеет первый порядок аппроксимации по $\tau$ и $h$ .

Замечание. Из процедуры интегрирования следует, что искомые значения $u^{n-1\!\not{\phantom{|}}\,\,2},\, p^{n-1\!\not{\phantom{|}}\,\,2}$ соответствуют среднеинтегральным величинам

$u^{n-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} u(x,\tau)\,dx\,,\qquad p^{n-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h} \int\limits_{x_{i-1}}^{x_{i}} p(x,\tau)\,dx\,.$

Приведем пример расчета по описанной схеме для частного случая рассматриваемой задачи (где $\rho_0=0,\!123;~ a_0=0,\!25$ ):

$\begin{gathered}\begin{aligned}& \frac{\partial u(x,t)}{\partial t}+ \frac{1}{\rho_0}\cdot \frac{\partial p(x,t)}{\partial x}=0,\\ & \frac{\partial p(x,t)}{\partial t}+\rho_0a_0^2\cdot \frac{\partial u(x,t)}{\partial x}=0,\end{aligned}\qquad 0<x<1,\quad 0<t<1,\\[4pt] \begin{aligned}& u(x,0)= \psi_1(x)=\sin x, & & p(x,0)=\psi_2(x)= \rho_0x\sin x, & & 0 \leqslant x \leqslant 1,\\ & p(0,t)= \varphi_1(t)=0,\!1t+ \rho_0a_0^2\sin t, & & p(1,t)= \varphi_2(t)= 0,\!1t+ \rho_0a_0^2\sin t, & & 0<t \leqslant1, \end{aligned} \end{gathered}$

Решение

А. Выберем число промежутков $n=2$ . Тогда находится шаг $h=\frac{1}{2}= 0,\!5$ и вводятся узлы $x_0=0;~ x_1=0,\!5;~ x_2=1$ . Зададим шаг $\tau= 0,\!5$ по времени из условия $\tau< \frac{h}{a_0}= \frac{0,\!5}{0,\!25}=2$ . Сформируем две ячейки "A" и "B" первого временного слоя $t_1=\tau= 0,\!5$ (рис. 8.27).

Б. Заменим функции $\psi_1(x)=\sin x,~ \psi_2(x)=\rho_0x\sin x$ кусочно-постоянными по формулам (8.81). Для вычисления интегралов используются квадратурные формулы трапеций:

$u_{i-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h}\cdot \frac{u_{i-1}+ u_{i}}{2}\cdot h=\frac{u_{i-1}+ u_{i}}{2}\,,\qquad p_{i-1\!\not{\phantom{|}}\,\,2}= \frac{1}{h}\cdot \frac{p_{i-1}+ p_{i}}{2}\cdot h=\frac{p_{i-1}+ p_{i}}{2}\,,\quad i=1;2,$

где $u_{i}=\sin x_{i},~ p_{i}=\rho_0x_{i}\sin x_{i}$ . Поэтому

$\begin{aligned}& u_{1\!\not{\phantom{|}}\,\,2}\cong \frac{0+0,\!479}{2}\cong 0,\!239; & & u_{3\!\not{\phantom{|}}\,\,2}\cong \frac{0,\!479+ 0,\!841}{2}\cong 0,\!660;\\[2pt] <br />& p_{1\!\not{\phantom{|}}\,\,2}\cong \frac{0+0,\!0294}{2}\cong 0,\!239; & & p_{3\!\not{\phantom{|}}\,\,2}\cong \frac{0,\!479+ 0,\!841}{2}\cong 0,\!660<br />\end{aligned}$

В. Найдем значения $P_0,\,P_2$ по формулам (8.82). Для вычисления интегралов, как и ранее, применяется формула трапеций:

$P_0\cong \frac{0+0,\!0537}{2}\cong 0,\!0269;\qquad P_2\cong \frac{0+0,\!0463}{2}\cong 0,\!0231.$

Рассчитаем значение $U_0$ на левой границе ячейки "А" по формуле (8.86): $U_0=\frac{P_0}{\rho_0a_0}+ u_{1\!\not{\phantom{|}}\,\,2}-\frac{p_{1\!\not{\phantom{|}}\,\,2} }{\rho_0a_0}\cong 0,\!635$ и значение $U_2$ на правой границе ячейки по формуле (8.87): $U_2=-\frac{P_2}{\rho_0a_0}+ u_{3\!\not{\phantom{|}}\,\,2}+ \frac{p_{3\!\not{\phantom{|}}\,\,2} }{\rho_0a_0}\cong 1,\!476$ .

Г. Вычислим значения $U_1,\,P_1$ на правой границе ячейки "А" (левой границе ячейки "В" ). По формулам (8.85) получаем

$U_1=\frac{u_{1\!\not{\phantom{|}}\,\,2}+ u_{3\!\not{\phantom{|}}\,\,2}}{2}-\frac{p_{3\!\not{\phantom{|}}\,\,2}-p_{1\!\not{\phantom{|}}\,\,2}}{2\rho_0a_0}\cong-0,\!392;\qquad \frac{p_{1\!\not{\phantom{|}}\,\,2}+ p_{3\!\not{\phantom{|}}\,\,2}}{2}-\rho_0a_0\, \frac{u_{3\!\not{\phantom{|}}\,\,2}-u_{1\!\not{\phantom{|}}\,\,2}}{2}\cong 0,\!0341.$

Д. Вычислим значения $u^{1\!\not{\phantom{|}}\,\,2},\, p^{1\!\not{\phantom{|}}\,\,2},\, u^{3\!\not{\phantom{|}}\,\,2},\, p^{3\!\not{\phantom{|}}\,\,2}$ на верхних границах ячеек:

$\begin{aligned}& u^{1\!\not{\phantom{|}}\,\,2}= u_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_1-P_0)\cong 0,\!1805; & & p^{1\!\not{\phantom{|}}\,\,2}= p_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0a_0^2(U_1-U_0)\cong 0,\!0226;\\[4pt] & u^{3\!\not{\phantom{|}}\,\,2}= u_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_2-P_1)\cong 0,\!749; & & p^{3\!\not{\phantom{|}}\,\,2}= p_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0a_0^2(U_2-U_1)\cong 0,\!0521. \end{aligned}$

Таким образом, первый слой полностью рассчитан.

Е. Для расчета второго слоя $(t_2=t_1+\tau)$ с шагом $\tau=0,\!5$ полученные значения $u^{1\!\not{\phantom{|}}\,\,2},\, p^{1\!\not{\phantom{|}}\,\,2},\, u^{3\!\not{\phantom{|}}\,\,2},\, p^{3\!\not{\phantom{|}}\,\,2}$ переобозначим, для чего верхние индексы заменяются нижними:

$u_{1\!\not{\phantom{|}}\,\,2}=0,\!1805;\quad p_{1\!\not{\phantom{|}}\,\,2}= 0,\!0226;\quad u_{3\!\not{\phantom{|}}\,\,2}=0,\!749;\quad p_{3\!\not{\phantom{|}}\,\,2}= 0,\!0521$

Сформируем две ячейки "С" и "D" второго временного слоя (см. рис. 8.27).

Ж. Найдем значения $P_0,\,P_2$ по формулам (8.82). Для вычисления интегралов, как и ранее, применяется формула трапеций:

$P_0\cong \frac{0,0537+ 0,1065}{2}\cong 0,\!0801;\qquad P_2\cong \frac{0,0463+ 0,00935}{2}\cong 0,\!0278.$

Рассчитаем значение $U_0$ на левой границе ячейки "С" по формуле (8.86):

$U_0= \frac{1}{\rho_0\,a_0}\cdot P_0+ u_{1\!\not{\phantom{|}}\,\,2}-\frac{1}{\rho_0\, a_0}\cdot p_{1\!\not{\phantom{|}}\,\,2}\cong 2,\!0504$

и значение $U_2$ на правой границе ячейки "D" по формуле (8.87):

$U_2=-\frac{1}{\rho_0\,a_0}\cdot P_2+ u_{3\!\not{\phantom{|}}\,\,2}+ \frac{1}{\rho_0\,a_0}\cdot p_{3\!\not{\phantom{|}}\,\,2}\cong 0,\!8651.$

3. Вычислим значения $U_1,\,P_1$ на правой границе ячейки "С" (левой границе ячейки "D"). По формулам (8.85) получаем

$U_1= \frac{u_{1\!\not{\phantom{|}}\,\,2}+ u_{3\!\not{\phantom{|}}\,\,2}}{2}-\frac{p_{3\!\not{\phantom{|}}\,\,2}-p_{1\!\not{\phantom{|}}\,\,2}}{2\,\rho_0\,a_0}\cong-0,\!0150;\qquad P_1=\frac{p_{1\!\not{\phantom{|}}\,\,2}+ p_{3\!\not{\phantom{|}}\,\,2}}{2}-\rho_0a_0\, \frac{u_{3\!\not{\phantom{|}}\,\,2}-u_{1\!\not{\phantom{|}}\,\,2}}{2}\cong 0,\!02861.$

Вычислим значения $u^{1\!\not{\phantom{|}}\,\,2},\, p^{1\!\not{\phantom{|}}\,\,2},\, u^{3\!\not{\phantom{|}}\,\,2},\, p^{3\!\not{\phantom{|}}\,\,2}$ на верхних границах ячеек "С" и "D":

$\begin{aligned}& u^{1\!\not{\phantom{|}}\,\,2}= u_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_1-P_0)\cong 0,\!599; & & p^{1\!\not{\phantom{|}}\,\,2}= p_{1\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0\,a_0^2(U_1-U_0)\cong 0,\!0385;\\[2pt] & u^{3\!\not{\phantom{|}}\,\,2}= u_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h\,\rho_0}(P_2-P_1)\cong 0,\!7556; & & p^{3\!\not{\phantom{|}}\,\,2}= p_{3\!\not{\phantom{|}}\,\,2}-\frac{\tau}{h}\,\rho_0\,a_0^2(U_2-U_1)\cong 0,\!04533. \end{aligned}$

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

Замечание. Здесь рассмотрен пример методики конструирования метода С.К. Годунова и ее реализации для относительно простой линейной задачи. В нелинейном случае метод существенно усложняется. Для более глубокого изучения проблемы авторы рекомендуют, где данный метод применяется для различных классов нелинейных задач и задач с тремя и четырьмя независимыми переменными.

Математический форум (помощь с решением задач, обсуждение вопросов по математике).

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

Список статей » Список форумов

Часовой пояс: UTC + 3 часа [ Летнее время ]