8. Уравнения в частных производных Уравнение, связывающее неизвестную функцию, независимые переменные и частные производные неизвестной функции, называется уравнением в частных производных. Порядок n старшей производной называется порядком дифференциального уравнения. Уравнение в частных производных называется линейным, если оно линейно относительно неизвестной функции и её частных производных. Приведем основные уравнения математической физики, которые являются линейными уравнениями в частных производных второго порядка. 1. Волновое уравнение (уравнение колебаний) (8.1) описывает различные виды волн звуковые, упругие, электромагнитные и другие колебательные процессы. Функция зависит от пространственных переменных x, y, z и времени t.

2. Процесс распространения тепла в однородном изотропном теле описывается уравнением теплопроводности: (8.2) Уравнение теплопроводности в общем виде записывается так: (8.3) Здесь температура в точке (x, y, z) в момент времени t, теплоёмкость в точке (x, y, z) в момент времени t, коэффициент теплопроводности в точке (x, y, z) в момент времени t, плотность источников тепла. 3. Установившееся тепловое состояние в однородном изотропном теле описывается уравнением Пуассона (8.4)

Установившееся тепловое состояние в однородном изотропном теле при отсутствии источников тепла внутри тела описывается уравнением Лапласа (8.5) В уравнениях (8.4) и (8.5) функция зависит только от пространственных переменных x, y, z. Согласно классификации уравнений второго порядка уравнение (8.1) относится к гиперболическому типу, уравнение (8.2) к параболическому типу, а уравнения (8.4), (8.5) к эллиптическому типу. В конкретной постановке задачи математической физики необходимо найти решение одного из уравнений (8.1) (8.5), удовлетворяющее дополнительным (начальным и граничным) условиям. Начальные условия задаются с уравнениями (8.1) (8.3) и обычно имеют вид (8.6) (8.7) При этом для уравнения (8.1) искомое решение должно удовлетворять в начальный момент времени t 0 обоим условиям (8.6) и (8.7), а для уравнения (8.2) или (8.3) одному из условий (8.6), (8.7). Граничные условия для уравнений (8.1) (8.3): искомое решение должно удовлетворять на границе тела (или среды) одному из условий

(8.8) (8.9) где производная по нормали к границе тела. Для уравнений Пуассона или Лапласа задаются только граничные условия (8.10) или (8.11) Задача для уравнения Лапласа (8.5) с граничным условием (8.10) называется задачей Дирихле, а с условием (8.11) задачей Неймана. Область, описывающая тело, может быть задана в трехмерном, двухмерном или одномерном пространствах. В первом случае функция зависит от трех пространственных переменных и времени t, во втором от двух пространственных переменных и времени t:, а в третьем случае от переменных x и t:.

8.1.Разностный метод для уравнения колебаний Рассмотрим задачу о малых колебаниях натянутой струны с распределенной по длине нагрузкой f(x, t) (см. рис. 8.1): (8.11) (8.12) (8.13) Рис Уравнение колебаний струны. Явная схема 0 a x u t T u(x, t) μ(x) μ(x) μ2(t) μ2(t) μ1(t) μ1(t)

Струна совершает плоские колебания, т.е. точки струны перемещаются параллельно плоскости t = 0. Функция u(x, t) выражает смещение точки x струны в момент времени t от прямолинейной формы. Начальные условия (8.12) означают следующее. Форма струны в начальный момент времени t = 0 выражается функцией μ(x). Скорость перемещения точки x струны в момент времени t = 0 равна значению функции μ 0 (x). Краевые условия (8.13) говорят о том, что левый конец струны с течением времени совершает смещение μ 1 (t), а правый конец смещение μ 2 (t). Если концы струны закреплены, то μ 1 (t) = μ 2 (t) = 0. Мы предполагаем, что начальные условия (8.12) и краевые условия (8.13) должны быть согласованы между собой в угловых точках, т.е. выполнены условия. На рис. 8.1 представлен случай, когда, Введем сеточную область (рис.8.2, a)). В прямоугольной области зададим точки: (8.14)

Рассмотрим уравнение (8.11) в точках,,, и заменим производные разностными формулами (8.15) (8.16) Обозначим через приближенные значения искомой функции в точках. Тогда из уравнения (8.11) получим разностное уравнение, которое аппроксимирует уравнение (8.11) с порядком O(h 2 + τ 2 ): (8.17) Рис t x x1x1 t1t1 t2t2 t3t3 a)a) b)b) x i, t k–1 x i, t k+1 x i+1, t k x i–1, t k x2x2 x3x3

На рис. 8.2 b) изображен шаблон «крест» разностного уравнения (8.17). Разностное уравнение (8.17) связывает значения неизвестной функции на трех слоях (k – 1, k, k + 1). На слое k = 0 заданы начальные условия (8.12), из которых следует, что Чтобы найти значения неизвестной функции на слое k = 1, используем условие для производной ut(x, 0) из (8.12). Для этого построим разложение в ряд Тейлора (8.18) Из уравнения (8.11), учитывая первое условие в (8.12), выразим вторую производную (8.19) Теперь, учитывая условие в (8.12), из (8.18), (8.19) выводим формулу для вычисления значений функции на первом слое: (8.20) учетом (8.13), окончательно получим для приближенных значений искомой функции на первом слое формулы (8.21)

Учитывая граничные условия (8.13) из (8.17) выводим формулы для вычисления значений на слоях : Мы получили явные формулы решения разностной задачи. Решение задачи называется устойчивым, если малым изменениям исходных данных отвечают малые изменения решения. Для выполнения условия устойчивости разностной схемы (8.17) необходимо и достаточно, чтобы выполнялось условие Куранта c < h. Алгоритм решения задачи для уравнения колебаний струны с помощью явной разностной схемы. 1. Построить сеточную область (8.21) выбирая шаги h и τ так, чтобы выполнялось условие устойчивости (условие Куранта) c < h.

2. Вычислить значения u i,k искомой функции для k = 0, 1: (8.22) (8.23) 3. Вычислить значения u i,k для : (8.24)

Уравнение колебаний струны. Неявная схема Построим неявную схему для уравнения колебаний струны (8.11): (8.25) Для устойчивости схемы (8.25) параметры c, h, τ, σ должны удовлетворять условию Если ¼ σ 1/2, то схема (8.25) безусловно устойчива. Шаблон схемы (8.25) изображен на рис. 8.3.

Рис. 8.3 Значения решения на нулевом и первом слоях вычисляют по формулам (8.22), (8.23). На каждом k-ом слое (k = 2, 3, …, M) решают методом прогонки систему уравнений относительно : (8.26) t x

Алгоритм решения неявной разностной схемы для уравнения колебаний струны. 0. Построить сеточную область, выбирая шаги h, τ и параметр σ так, чтобы выполнялось условие устойчивости Если ¼ σ 1/2, то можно выбирать h, τ произвольно. 1. Вычислить значения ui,k искомой функции для k = 0, 1: (8.27) (8.28) 2. Значения u i,k+1 для каждого k > 0 находим методом прогонки, последовательными вычислениями в несколько этапов Вычислим правые части (8.26): (8.29)

2.2. Вычислим прогоночные коэффициенты: (8.30) (8.31) (8.32) 2.3. Вычислим решение u i,k+1 : (8.33) (8.34)

Разностный метод для уравнения колебаний мембраны Рассмотрим задачу для уравнения колебаний однородной прямоугольной мембраны (8.37) (8.38) (8.39) Введем сеточную область Обозначим. Заменяя производные разностными формулами, для уравнения (8.37) получим разностное уравнение с порядком аппроксимации :

(8.40) Алгоритм решения явной разностной схемы для уравнения колебаний мембраны. 1. Решение на нулевом слое k = 0 получим из начального условия (8.38): (8.41) Чтобы найти решение на первом слое (k = 1) используем разложение в ряд Тейлора по переменной t в окрестности t = 0: 2. Учитывая начальные и краевые условия, получим приближенные значения решения на первом слое: (8.42)

3. Для вычисления значений на следующих слоях (k > 0) имеем из (8.40): (8.43) Приведем условие устойчивости явной схемы (8.41) (8.43) (8.44)

8.2.Разностный метод для уравнения теплопроводности Рассмотрим краевую задачу линейной теплопроводности в однородной среде: (8.45) (8.46) (8.47) Краевые условия (8.47) называются условиями первого рода и означают, что на концах отрезка [0, a] задана зависимость температуры u от времени t. Краевые условия второго рода (8.48) означают, что на концах отрезка [0, a] заданы тепловые потоки. Если на границах имеется линейный теплообмен с окружающей средой, то задаются краевые условия третьего рода (8.49)

Одномерное уравнение. Неявная схема Введем сеточную область: Построим для задачи (8.45) (8.47) неявную схему с шаблоном, который изображен на рис. 8.4: (8.50) Если σ = 0, то (8.50) представляет явную схему: (8.51)

Порядок аппроксимации схемы (8.50) при σ ½ равен O(τ + h2). При σ = ½ схема (8.50) называется схемой с полусуммой или симметричной и имеет порядок сходимости O(τ2 + h2). Рис. 8.4 Из краевых условий (8.47) находим (8.52) Преобразуем разностную схему (8.50), (8.52) к виду (8.53) Для простоты записи в (8.53) введено обозначение: (8.54) t x

Разностная схема (8.53) решается методом прогонки. Положим σ = ½. Тогда (8.53), (8.54) переписываются в следующем виде: (8.55) (8.54) Построим алгоритм решения задачи (8.45) (8.47) с помощью неявной схемы (8.50) при σ = ½. 1. Из начального условия (8.46) находим значения неизвестной функции на нулевом слое: (8.56) 2. Значения функции на следующих слоях находим методом прогонки. Для k = 0, 1, …, M – 1 выполняем следующие вычисления:

2.1. Вычисляем правые части (8.55) по формулам: (8.57) (8.58) 2.2. Вычисляем коэффициенты прогонки: (8.59) (8.60) (8.61) 2.3. Вычислим решение u i,k +1 : (8.62) (8.63)

Двумерное уравнение. Неявная схема Рассмотрим первую краевую задачу для двумерного уравнения теплопроводности с постоянным коэффициентом в прямоугольной области: (8.65) (8.66) (8.67) Введем сеточную область

Обобщим неявную разностную схему (8.50) на двумерный случай (8.68) Шаблон схемы (8.68) изображен на рис Запишем для сеточной функции условия, вытекающие из (8.66), (8.67). Из (8.66) получим значения на нулевом слое: (8.69) Краевые условия (8.67) записываются в виде: (8.70)

Рис. 8.5 Разностная схема (8.68) (8.70) имеет порядок аппроксимации при σ ½, и при σ = ½. Схема устойчива при. (8.71) Очевидно, при σ ½ схема абсолютно устойчива. Для решения системы (8.68) (8.70) можно применить метод Гаусса на каждом слое k = 1, 2, …, M. t x y a b

Двумерное уравнение. Экономичные схемы На решение системы (8.68) (8.70) методом Гаусса потребуется порядка N 3p операций [9], где p размерность задачи, а N число интервалов по каждой пространственной переменной. В двумерном случае имеем N 6 операций, а в случае трехмерной задачи В двумерном случае имеем N 6 операций, а в случае трехмерной задачи число операций составит порядка N 9. Для многомерного параболического уравнения построены так называемые экономичные схемы, которые требуют для вычисления значений функции на одном слое порядка N 2 операций в двумерном случае и N 3 операций в трехмерном случае и являются абсолютно устойчивыми. Большинство расчетов для многомерного параболического уравнения проводится по таким схемам. Рассмотрим одну из таких схем, схему переменных направлений (продольно-поперечную схему). Выберем шаблон, представленный на рис. 8.6, содержащий промежуточный слой, и построим следующую схему:

(8.72) (8.73) Порядок аппроксимации схемы (8.72), (8.73) на целых слоях равен а на полуслоях составляет. Рис. 8.6 Начальные условия (8.66) дают значения сеточной функции на нулевом слое: (8.74) t x y tk tk t k + τ

Значения на промежуточном слое определяются методом прогонки из системы уравнений (8.75) Значения сеточной функции при находим из граничных условий (8.67) Однако здесь мы не будем использовать очевидные формулы так как они имеют порядок аппроксимации. Мы воспользуемся формулами, полученными следующим способом. Вычтем (8.73) из (8.72): (8.76)

Отсюда при i = 0, учитывая (8.67), получим (8.77) Аналогично при i = N x получим (8.78) Эти формулы имеют порядок.

Значения сеточной функции на следующем слое определяются методом прогонки из системы уравнений (8.79) (8.80) Алгоритм метода переменных направлений. 0. При k = 0 (8.81) 1. Для каждого значения k = 0, 1, …, M – 1 выполняем следующие два этапа вычислений: I –й этап. Выполняем для следующие 3 шага: 1.1. Вычисляем для : (8.82)

1.2. Вычисляем прогоночные коэффициенты: (8.83) (8.84) (8.85) (8.86)

(8.87) (8.88) 1.3. Вычислим решение : (8.89) (8.90) (8.91) (8.92) II-й этап. Выполняем для следующие 3 шага: 2.1. Вычисляем для : (8.93)

2.2. Вычисляем прогоночные коэффициенты: (8.94) (8.95) (8.96) 2.3. Вычислим решение : (8.97) (8.98) (8.99) (8.100) (8.101) (8.102)

8.3 Разностный метод для эллиптического уравнения. Рассмотрим задачу Дирихле для уравнения Пуассона в прямоугольной области (8.103) (8.104) Введем сеточную область Рассмотрим уравнение (8.103) в узлах сетки и заменим производные разностными формулами. Тогда с учетом краевых условий получим Задача Дирихле для уравнения Пуассона

систему уравнений для определения приближенного решения задачи (8.103), (8.104): (8.105) (8.106) Система (8.105), (8.106) содержит неизвестных значений и столько же уравнений. Матрица этой системы имеет блочно-трехдиагональную структуру. Перепишем систему (8.105), (8.106) в виде системы векторных уравнений Решение разностной краевой задачи для уравнения Пуассона методом Гаусса

(8.107) Каждое векторное уравнение в (8.107) содержит линейных уравнений, A и С матрицы размерности, Uj, Fj векторы-столбцы с числом элементов : (8.108) (8.109)

(8.110) (8.111)

(8.112) Система уравнений (8.107) может решаться следующим методом матричной прогонки: 1. Вычисляем векторы F j по формулам (8.110) (8.112). 2. Вычислим матрицы A и С по формулам (8.108). 3. Вычисляем прогоночные коэффициенты-матрицы, : (8.113) 4. Вычисляем решение системы: (8.114)

Метод матричной прогонки легко реализуется в виде программы для компьютера. Необходимость вычисления на каждом шаге в (8.113) обратной матрицы можно считать недостатком метода, так как обращение матрицы требует выполнения большого числа операций. Этого недостатка можно избежать, применяя метод исключений Гаусса с учетом разреженности матриц A и С в системе (8.107). Приведем общее описание алгоритма метода Гаусса: 1. Прямой ход. С помощью обычных исключений Гаусса, преобразуем первое векторное уравнение системы (8.107) к виду, где C 1 верхняя треугольная матрица с единицами на главной диагонали. С помощью соответствующих уравнений из векторного уравнения, исключаем слагаемое из второго векторного уравнения. Затем, полученное второе векторное уравнение преобразуется к виду, C 2 верхняя треугольная матрица с единицами на главной диагонали. Продолжая этот процесс, приводим систему (8.107) к виду (8.115) где C j верхние треугольные матрицы с единицами на главной диагонали.

2. Обратный ход. Из последнего векторного уравнения в системе (8.115) находим вектор, а затем последовательно вычисляем для с помощью первого уравнения в (8.115) Решение разностной краевой задачи для уравнения Пуассона методом переменных направлений Задача (8.103), (8.104) является стационарной задачей, её решение не зависит от времени. Рассмотрим эволюционную задачу для параболического уравнения с такими же краевыми условиями (8.104) и произвольно заданными начальными условиями. (8.116) (8.117) (8.118)

Приведем без доказательства следующий факт : Если начальные и краевые условия задач (8.103), (8.104) и (8.116) (8.118) таковы, что их решения имеют в прямоугольнике непрерывные производные, равномерно ограниченные по t, то при решение эволюционной задачи (8.116) (8.118) равномерно сходится по t к решению стационарной задачи (8.103), (8.104). Физический смысл этого факта заключается в следующем. Задача (8.116) (8.118) описывает изменение с течением времени температуры в точках прямоугольной области при заданных температурах в начальный момент времени и постоянных значениях температуры на границе прямоугольника. С течением времени процесс изменения температуры выходит на стационарный режим, т.е. устанавливается решение перестает зависеть от времени. Задача (8.103), (8.104) как раз и соответствует установившемуся процессу. Этот факт позволяет построить эффективный итерационный метод решения задачи (8.103), (8.104), который заключается в том, чтобы найти решение задачи (8.116) (8.118) при больших значениях t, т.е. найти решение, не зависящее от t. Такой способ решения задачи для эллиптического уравнения называется счетом на установление.

Таким образом, для решения задачи (8.103), (8.104), можно применить метод переменных направлений к задаче (8.116) (8.118), причем надо вычислять до такого значения k, чтобы выполнилось неравенство (8.119) где ε заданная точность. При этом в алгоритме (8.81) (8.102) следует учесть независимость граничных условий от времени и выбрать произвольные начальные условия.