Уравнения математической физики

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

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

, (7.1)

где F – заданная функция, – искомая функция.

Далее будут рассматриваться уравнения более простого вида, а именно линейные уравнения первого и второго порядков с постоянными коэффициентами. Такие уравнения имеют вид:

. (7.2)

В общем же случае коэффициенты a, b, c, d, e, f также как и правая часть g, могут зависеть от переменных x, y и искомой функции u.

Существуют различные виды уравнений в зависимости от соотношения между коэффициентами. Рассмотрим некоторые из них. При a = b = c= = f = 0, , получается уравнение первого порядка вида

,

называемое уравнением переноса. На практике в этом уравнении одной из переменных может быть время t. Тогда его называют также эволюционным уравнением.

Если хотя бы один из коэффициентов a, b, c отличен от нуля, то (7.2) является уравнением второго порядка. В зависимости от знака выражения оно может принадлежать к одному из трех типов: гиперболическому ( ), параболическому ( ) или эллиптическому ( ).

Приведем примеры уравнений второго порядка с частными производными разных типов:

– волновое уравнение (гиперболическое)

;

– уравнение теплопроводности или диффузии (параболическое)

, ;

– уравнение Пуассона (эллиптическое)

.

Если правая часть последнего уравнения равна нулю , то оно называется уравнением Лапласа.

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

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

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

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

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

Решение дифференциальных уравнений рассматривается, как правило, для систем, ограниченных в пространстве и в течение конечного промежутка времени. Следовательно, решение уравнения ищется в ограниченной, в данном случае прямоугольной, области G:

. (7.3)

Для получения частного решения параболического уравнения (уравнения теплопроводности) необходимо задать одно начальное условие и два граничных. Начальное условие может быть задано в виде

, , (7.4)

где – заданная функция, определяющая значение функции в начальный момент времени .

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

, , , (7.5)

где и – заданные функции.

Например, если – температура тонкого однородного стержня длинной l, то условия (7.4) задает начальное значение температуры каждой точки стрежня при , а условия (7.5) означают, что температура в начале и конце стержня задана в каждый момент времени, в частности, поддерживается постоянной и .

Геометрической интерпретацией искомого решения является поверхность в пространстве , которая проецируется на область G плоскости , причем в соответствии с (7.4) и (7.5) заданы три кривые, которые представляют собой края этой поверхности (рис. 7.1). Внутри области G значения функции неизвестны, поэтому неизвестна форма поверхности и ее требуется найти путем решения начально-краевой нестационарной задачи.

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

, , , (7.8)

где и – заданные функции, определяющие значение функции и ее первой производной по t в начальный момент времени . Граничные условия задаются аналогично (7.5). Геометрически решению данной начально-краевой задачи соответствует поверхность в пространстве (рис. 7.1).

§ 2. Разностные методы решения уравнений
с частными производными

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

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

В случае прямоугольной области G расчетная сетка строится следующим образом (рис. 7.2). Отрезок на осиx разбивается на n равных частей длины , при этом получается узловых точек , . Аналогично отрезок на оси y разбивается на m равных частей длины , при этом получается узловых точек , . Проводя через эти точки прямые, параллельные осям координат, получаем сетку, разбивающую область G на элементарные прямоугольные ячейки. Любой узел сетки, номер которого , определяется координатами . Решение задачи ищется именно в этих точках. Поскольку все ячейки построенной сетки одинаковы, такую сетку называют равномерной. Узлы сетки, лежащие на границе области G, называютсяграничными узлами; все остальные узлы – внутренними.

Аналогично вводятся сетки и для многомерных областей, содержащих более двух измерений.

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

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

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

, , , ,

(7.9)

, , ,

где – начальное распределение температуры U вдоль стержня (при ); , – значения температуры на концах стержня в любой момент времени t. Коэффициент a определяется теплопроводностью k и удельной теплоемкостью l материала стержня ( ). Заметим, что начальные и граничные условия должны быть согласованы, т.е. , .

Введем равномерную сетку с помощью координатных линий ( ), ( ); h и t – шаги сетки по направлениям x и t, соответственно. Значения искомой функции в узлах сетки обозначим . Эти значения заменим соответствующими значениями сеточной функции , которые удовлетворяют уравнениям, образующим разностную схему.

Заменим в уравнении (7.9) частные производные искомой функции приближенными разностными отношениями:

,

Подстановка этих соотношений в (7.9) дает систему разностных уравнений для внутренних узлов сетки

,

(7.10)

,

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

, ,

(7.11)

, ,

Соотношения (7.10) в совокупности с дополнительными условиями (7.11) образует замкнутую систему уравнений называемую разностной схемой.

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

Совокупность узлов при фиксированном значении j называется слоем (в данном случае временным слоем).

Схема (7.10) позволяет последовательно находить значения ( ) на -м слое через соответствующие значения на j-том слое:

. (7.12)

Такие схемы называются явными.

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

то вместо (7.10) получим разностную схему

,

(7.13)

, .

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

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

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

.

В качестве глобальной погрешности (т.е. погрешности решения на всей области G) можно принять максимальное по модулю значение на сетке

.

Разностная схема называется сходящейся, если при сгущении узлов сетки (т.е. при и ) значение погрешности стремится к нулю:

.

Если при этом ), где M – некоторая положительная константа, то разностная схема имеет p-й порядок точности по h и q-й порядок точности поt. Говорят также, что она сходится со скоростью .

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

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

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

. (7.14)

При этом для погрешности справедлива оценка

.

То есть явная схема (7.10) имеет второй порядок точности по h и первый порядок точности по t.

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

В тоже время доказано, что неявная схема (7.13) является безусловно устойчивой.

Условие устойчивости (7.14) явной схемы накладывает существенные ограничения на выбор шага t по времени. Так, например, если пространственный шаг hимеет порядок , а коэффициент a – порядок единицы, то t ограничено величиной порядка . Следовательно, может потребоваться очень большое количество шагов по времени для того, чтобы получить решение во всей расчетной области G. Применение неявной схемы, позволяет проводить расчет с существенно бóльшим значением шага по времени и, следовательно, решение может быть получено за меньшее количество шагов. Однако погрешность получаемого решения, также как и для явной схемы, определяется величиной , которая содержит t в первой степени и поэтому не позволяет получать решения высокой точности при слишком больших значениях шага по времени.

Волновое уравнение. Уравнение гиперболического типа, называемое также волновым уравнением, возникает при описании различных колебательных процессов. Поскольку колебания представляют собой нестационарный процесс, то одной из независимых переменных является время t. Другими независимыми переменными выступают пространственные координаты x, y и z. В зависимости от их количества различают одномерное, двумерное и трехмерное волновые уравнения.

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

В качестве примера рассмотрим одномерное волновое уравнение

, (7.15)

которое описывает свободные поперечные колебания упругой натянутой струны. Искомая функция описывает положение (форму) струны в момент времениt. Коэффициент a определяется натяжением струны T и ее линейной плотностью r: . Колебания предполагаются малыми, т.е. их амплитуда мала по сравнению с длинной струны.

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

, , (7.16)

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

Если концы струны жестко закреплены, то их смещения равны нулю и граничные условия имею вид:

, . (7.17)

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

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

, (7.18)

. (7.19)

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

, . (7.20)

Соответствующий этой схеме шаблон показан на рис. 7.4.

Из (7.20) можно найти явное выражение для значений сеточной функции на (j+1)-м слое:

, . (7.21)

Полученная явная схема является трехслойной. Это означает, что для определения неизвестных значений на (j+1)-м слое нужно знать решения на j-м и (j-1)-м слоях. Следовательно, для того чтобы начать счет по формулам (7.21) необходимо знать решения на нулевом и первом временных слоях. Их можно найти исходя из начальных условий (7.16).

Так для нулевого слоя имеем

, (7.22)

Для получения решения на первой слое привлечем второе начальное условие (7.16). Заменим первую производную приближенным разностным отношением

, (7.23)

Из этого соотношения можно выразить значения искомой функции на первом временном слое:

, . (7.24)

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

, , . (7.25)

Полученная расчетная схема (7.21) является условно устойчивой. Необходимое и достаточное условие устойчивости имеет вид

. (7.26)

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

. (7.27)

Первую производную по t можно заменить функцией . Выражение для второй производной можно найти использую исходное дифференциальное уравнение (7.15) и первое начальное условие (7.16):

. (7.28)

С учетом (7.27) и (7.28) получим приближенное соотношение для искомой функции на первом временном слое

, . (7.29)

Разностная схема (7.20) с учетом (7.29) обладает общей погрешностью порядка .

Помимо колебаний струны уравнение (7.15) описывает также распространение плоских волн в упругой среде, например в атмосфере или слое жидкости. Параметр a в этом случае имеет смысл скорости распространения волн, которая может достигать очень больших значений (нескольких километров в секунду в задачах гидроакустики). При этом условие устойчивости (7.26) накладывает сильное ограничение на выбор величины шага по времени t, что существенно увеличивает необходимое компьютерное время и затраты на проведение расчетов. Чтобы избежать ограничения на выбор шага t, связанного с необходимостью выполнения условия (7.26), можно использовать неявную разностную схему, которая строится следующим образом.

Заменим производные, входящие в исходное дифференциальное уравнение (7.15) в узле ( ) разностными отношениями, построенными по пятиточечному шаблону, показанному на рис. 7.5.

,(7.30)

(7.31)

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

, . (7.32)

Систему (7.32) необходимо дополнить соотношениями (7.22), (7.24) и (7.25), которые следуют из начальных и граничных условий. Получившуюся систему алгебраических уравнений можно решить, последовательно определяя значения сеточной функции на слоях , методом прогонки (см. гл. 3). Действительно, при соотношение (7.32) принимает вид

,

где значения и известны из начальных условий, поэтому получаем систему с трехдиагональной матрицей относительно неизвестных , которая может быть решена методом прогонки. Значения в крайних точках и определяются из граничных условий (7.25). После нахождения приближенного решения на втором слое, аналогичной находится решение на третьем слое и т.д. Применение неявной разностной схемы (7.32) не требует выполнения какого-либо соотношения между шагами h и t, которое бы обеспечивало устойчивость приближенного решения, что позволяет производить вычисления с большим шагом по времени, чем при использовании явной разностной схемы (7.20).

Уравнение Лапласа. Задача Дирихле. Многие стационарные физические задачи (т.е. такие задачи, в которых рассматриваются явления, независящие от времени) сводятся к решению уравнения Пуассона вида

. (7.33)

Если правая часть равная нулю , то уравнение (7.33) называется уравнением Лапласа. Далее для простоты будем рассматривать двумерное уравнение Лапласа

. (7.34)

Решение этого уравнения будем искать на ограниченной области G изменения независимых переменных x и y. Границей области является замкнутая линия L. В простейшем случае область может иметь прямоугольную форму. Для полной формулировки краевой задачи кроме самого уравнения (7.34) нужно задать граничное условие на границе L. Запишем его в следующем виде

. (7.35)

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

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

. (7.36)

Из (7.36) можно получить систему линейных алгебраических уравнений относительно неизвестных значений сеточной функции в виде

,

, . (7.37)

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

, , ;

, , ;

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

. (7.38)

На основе этого уравнения строят итерационную процедуру:

, (7.39)

Начальное приближение выбирается произвольным образом (иногда полагают для всех i, j внутри области G).

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

.

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

Поскольку решение уравнения Лапласа не зависит от времени то формально в него можно добавить нулевое слагаемое . Тогда уравнение (7.34) примет следующий вид

. (7.40)

Полученное уравнение представляет собой двумерное уравнение теплопроводности. Разностные схемы для этого уравнения строятся точно также, как это было сделано нами для одномерного уравнения теплопроводности. Это уравнение остается дополнить начальным условием, например таким:

. (7.41)

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