Численные методы решения дифференциальных уравнений в частных производных. 2 страница
· для отрезка [a, b]:
(27)
Формула Симпсона-Ньютона используется преимущественно в расчетах на компьютерах, так как вычисления по ней программируются так же легко, как и расчет по формуле трапеций, а погрешность вычисления существенно меньше (формула трапеций имеет третий порядок точности, а формула Симпсона-Ньютона – пятый).
3.5. Численное решение обыкновенных дифференциальных уравнений
Обыкновенные дифференциальные уравнения (ОДУ) и системы таких уравнений широко встречаются в математическом моделировании. ОДУ связывает между собой значения независимой переменной (аргумента), искомой величины – функции этого аргумента и производных этой функции. Наибольший порядок производной определяет и порядок уравнения.
ОДУ первого порядка широко встречаются в ММ динамических дискретных систем, ОДУ второго и высших порядков могут возникать и в моделях континуальных стационарных систем, если характеристики системы непрерывно меняются вдоль одной пространственной координаты.
При интегрировании ОДУ появляется ряд неопределенных констант интегрирования, число которых совпадает с порядком уравнения. Для определения этих констант служат начальные и граничные условия.
Начальные условия – это задание определенных связей для искомой функции при начальном значении аргумента (в начальный момент времени, в начальной точке какой-то линейно протяженной пространственной области). Задача отыскания решения системы ОДУ первого порядка с заданными начальными условиями называется в математике задачей Коши.
Граничные условия задаются при двух значениях аргумента. Обычно это условия на границах пространственной области.
а) тело падает с некоторой высоты. Необходимо определить зависимость скорости и положения тела над землей от времени.
На тело действуют две силы: сила тяжести, равная mg, где m – масса тела, g– ускорение свободного падения, и сила сопротивления воздуха, имеющая две составляющих – силу трения, пропорциональную скорости тела, а также вихревое сопротивление, пропорциональное квадрату скорости. Не вдаваясь в физические подробности, запишем общее сопротивление так: Fсопр=-(av+bv 2 ).Теперь уравнение второго закона Ньютона применительно к телу в проекции на ось Oy запишется (для удобства направим ось Оу вниз):
. (28)
. (29)
Таким образом, движение тела описывается уравнением (штрих, как обычно в математике, обозначает производную):
. (30)
Оно имеет второй порядок и требует постановки двух начальных условий. Если отсчитывать координату у от начального положения тела и считать, что его движение начинается из состояния покоя, условия примут вид:
. (31)
Задачу можно описать также системой двух уравнений (28), (29), каждое из которых имеет первый порядок. Таким образом, уравнение высокого порядка можно преобразовать в систему уравнений пониженного порядка, принимая производные за новые переменные аналогично (29). Начальные условия в этом случае будут выглядеть:
. (32)
Рассмотренная задача дает пример математической модели физического процесса, описываемого задачей Коши.
б) для улучшения охлаждения к корпусу теплового двигателя температурой Т0 приварены ребра – тонкие металлические пластины. Необходимо определить распределение температуры вдоль ребра и количество тепла, отводимое с его помощью в окружающую среду.
Распространение тепла в толще ребра описывается так называемым законом Фурье:
, (33)
где Q – количество тепла, переносимое за 1 сек через сечение ребра площадью F=bH; l —так называемый коэффициент теплопроводности, Т – температура. Количество тепла, отводимое в окружающую среду с участка поверхности ребра длиной dx, равно (Тос – температура окружающей среды):
.
Первое слагаемое в квадратных скобках выражает охлаждение поверхности ребра омывающим его воздухом, второе – отвод тепла путем теплового излучения. Исходя из теплового баланса слоя ребра, ограниченного сечениями, соответствующими значениям координат x и x+dx, можем записать:
,
В левой части равенства стоит выражение производной
(так как наличие дифференциала dxуказывает на предельный переход при х®0, ср. (22)). С учетом (33) получаем уравнение:
. (34)
Это уравнение второго порядка и для него надо поставить два условия. Одно из них очевидно: T=T0при х=0. Второе условие в теории теплообмена ставится по-разному. Простейший вариант (пригодный для длинных ребер) – считать, что температура конца ребра совпадает с температурой окружающей среды: Т=Тоспри х=L. Таким образом, в данной задаче для дифференциального уравнения (34) поставлены граничные условия.
Рассмотрим некоторые методы решения ОДУ. Начнем с методов решения задачи Коши, то есть, численного решения дифференциальных уравнений вида
(35)
и систем уравнений вида
(36)
с начальными условиями
или 
(37)
Используя формулу производной сеточной функции (23), можно записать:
Это и есть расчетная формула метода Эйлера. Для решения отрезок оси Ох,на котором надо найти решение уравнения, разбивают на малые шаги h (создают сетку), затем в правую часть подставляют значения 
и 
. Это позволяет приближенно определить значение искомой функции y(x) в точке 
. Затем найденные значения 
и 
опять подставляют в правую часть формулы (6), вычисляют 
и т.д. Если обозначить 
, то формулу (6) можно записать в окончательном виде
(38)
б) метод Рунге-Кутта. Этот метод основан на использовании представлений производных, имеющих более высокий порядок точности. Вычисление решения в очередном узле проводится по следующим формулам (приводятся без вывода):
(39)
(40)
Оба метода достаточно легко обобщаются на случай систем ОДУ – в этом случае функции в правых частях формул оказываются функциями нескольких переменных, но их значения все равно вычисляются для уже рассчитанных узлов сетки.
Метод Эйлера реализуется существенно легче, чем метод Рунге-Кутта, поэтому его можно применить даже при расчете вручную, однако его точность крайне низка – он имеет первый порядок точности, тогда как метод Рунге-Кутта – четвертый. Поэтому в компьютерных расчетах используется метод Рунге-Кутта или другие методы высокого порядка точности.
Для решения краевых задач применяются разные методы. Рассмотрим самый простой из них – метод пристрелки или стрельбы. Он основан на преобразовании краевой задачи в задачу Коши.
Рассмотрим для примера задачу (34), которую запишем в виде:
(41)
Вместо последнего граничного условия при х=L зададим дополнительное условие при х=0: q = q0. Тогда мы получим задачу Коши. Величину q0 можно выбрать произвольно, например, принять в первом приближении равной нулю.
Решив задачу Коши одним из вышеописанных методов, мы определим значение неизвестной величины Т на конце ребра, при х=L. Скорее всего, оно будет отличаться от значения Тос. Сущность метода пристрелки заключается в том, что значение дополнительного начального условия подбирается таким образом, чтобы после решения задачи Коши оказалось выполненным граничное условие на конце рассматриваемого отрезка (в данном случае при x=L). Это действительно напоминает пристрелку при артиллерийской стрельбе, где дополнительному начальному условию можно сопоставить угол наклона ствола орудия, который меняется до тех пор, пока не будет выполнено граничное условие на другом конце – совпадение траектории полета снаряда с точкой местонахождения цели.
Фактически получается, что значение граничного условия на конце отрезка является функцией величины дополнительного начального условия: T½x=L = f(q0) . Тогда условие совпадения конечного условия с заданной величиной Тос запишется в виде некоторого нелинейного уравнения:
,
где функция определяется путем решения задачи Коши. Подбор граничного условия можно, таким образом, проводить с помощью описанных ранее методов решения нелинейных алгебраических уравнений, например, метода половинного деления (что, опять-таки, имеет прозрачную «артиллерийскую» аналогию – взятие цели в «вилку»).
3.6. Моделирование процессов, приводящих
к дифференциальным уравнениям в частных производных
При анализе химико-технологических процессов, нагрева и охлаждения тел и т.п., часто приходится иметь дело с величинами, значения которых меняются от одной точки тела к другой непрерывно. Так, если поместить головку сыра для хранения в холодильник, ее температура будет меняться постепенно. Сначала остынет поверхностный слой, тогда как центральная часть головки будет оставаться теплой. Постепенно охлажденная зона будет занимать все более толстый слой, но в центре головки температура долго будет оставаться более высокой, чем у поверхности. Изменение температуры внутри головки можно представить графиком Т(х), где х – координата, отсчитываемая вдоль диаметра головки. Очевидно, температура является функцией расстояния от центра головки, а также зависит от времени t:
Эта функция непрерывна, т.к. температура в сплошной среде не может измениться скачкообразно. Если вместо симметричной головки взять тело иной формы, температура в нем окажется функцией трех координат и времени:
(42)
Такое непрерывное распределение температуры по точкам в сплошной среде называется полем температуры. Поля могут быть образованы и другими величинами – давлением, концентрацией растворённого вещества и т.п.
Бывают векторные поля; например, поле скоростей в жидкости или газе. Если в жидкости движется какое-то тело (например, весло в воде), оно увлекает за собой часть жидкости; в то же время, жидкость вдали от тела остается неподвижной. Вектор скорости изменяется от точки к точке, постепенно уменьшаясь до нуля, так что
.
Такое векторное поле можно представить, как комбинацию трех скалярных полей проекций скорости на оси координат.
3.7. Дифференциальное уравнение теплопроводности
Рассмотрим процесс теплообмена в твердом теле. Количество тепла, передаваемое за некоторое бесконечно малое время dtчерез бесконечно малый элемент dF поверхности, нормальной к координатной оси Ox, пропорционально dtи dF, а также зависит от скорости изменения температуры в направлении оси Ox– чем быстрее изменяется температура, тем больше передается тепла. В большинстве случаев связь количества передаваемого тепла со скоростью изменения температуры вдоль Ox принимается линейной:
. (43)
 |
Коэффициент пропорциональности l называется коэффициентом теплопроводности; знак «минус» указывает на известный из физики факт – тепло передается в сторону убывания температуры.
называется плотностью потока тепла в направлении координаты x (численно величина qxсоответствует количеству тепла, передаваемому за 1 сек через 1 м 2 поверхности). Аналогично можно записать:
Величины qx, qy, qz можно рассматривать, как компоненты вектора
,
называемого вектором плотности потока тепла. Он определяет направление переноса тепла через данную точку пространства и количество тепла, переносимого за 1 сек в расчете на 1 м 2 площади поверхности, нормальной к этому направлению. Производные температуры по координатам тоже являются компонентами вектора, называемого градиентом температуры:
.
Для определения градиента температуры необходимо знать вид функции (42), поэтому он является характеристикой температурного поля. Аналогично вводятся градиенты других величин – давления, концентрации и т.п. (существуют и градиенты векторных величин, например, скорости, но их рассмотрение требует применения сложного математического аппарата).
В векторной форме закон теплопроводности записывается:
.
Это соотношение называется законом Фурье. Аналогичный вид имеет так называемый закон Фика, описывающий диффузию какого-либо вещества в сплошной среде:
,
где вектор
— вектор плотности потока массы i—компонента смеси веществ (модуль этого вектора показывает, сколько данного вещества проходит за 1 сек через малый элемент поверхности, нормальной к направлению переноса, в расчете на 1 м 2 площади поверхности), grad ci– градиент концентрации данного вещества, r— плотность среды, Di – коэффициент диффузии данного вещества.
В силу аналогии законов Фурье и Фика многие закономерности и расчетные формулы теплообмена оказываются справедливыми и для массообмена (диффузии).
Рассмотрим бесконечно малый элемент сплошной среды в виде параллелепипеда с ребрами dx, dy, dz. Через левую грань объема в него поступает количество теплоты dQx, описываемое формулой (43). Одновременно через правую грань из него уходит теплота
Если температура меняется вдоль оси Ох, значения производных 
на противоположных гранях неодинаковы. В итоге внутри рассматриваемого элемента остается тепло
Если dx– бесконечно малая величина (т.е., в конце вывода должен быть совершен предельный переход при 
), выражение в квадратных скобках является второй производной температуры по х (со знаком «минус»); тогда
.
Аналогично можно рассмотреть две другие пары противоположных граней. Тогда общее количество тепла, накопленное в объеме за время dt, равно
.
Накопление тепла вызывает нагрев вещества в объёме:
,
где 
— частный дифференциал температуры.
Приравняв два выражения для d 2 Q, получаем
(44)
— дифференциальное уравнение теплопроводности. Оно получено в предположении постоянства l и отсутствия внутреннего тепловыделения в сплошной среде (такое тепловыделение может быть обусловлено, например, протеканием химической реакции в веществе или другими физико-химическими процессами). Более общие формы уравнения теплопроводности рассматриваются в литературе по тепломассообмену.
Символом 
обозначен так называемый лапласиан температуры (сумма вторых частных производных по координатам). Лапласианы различных величин часто встречаются в дифференциальных уравнениях в частных производных.
Аналогичным образом (применяя законы физики к бесконечно малым объемам) можно получить и другие дифференциальные уравнения, описывающие процессы в сплошных средах. Примеры таких уравнений:
;
— волновое уравнение (уравнение распространения волн в сплошной среде):
,
где u – величина, волновое колебание которой рассматривается (давление, температура, электрический потенциал и др.), v– скорость распространения волны;
, (45)
описывающее стационарные распределения различных величин (температуры, электростатического потенциала и др.) в сплошной среде.
3.8. Краевые задачи для уравнений в частных производных
Уравнение теплопроводности применимо к анализу распространения тепла в любой системе – от куска сыра до стенки атомного реактора. Поэтому при рассмотрении конкретных задач к уравнению надо добавить описание особенностей задачи, отличающих ее постановку от других задач теплопроводности.
Эти особенности называют условиями однозначности. Сюда входит описание формы и указание размеров тела, указание значений его физических свойств.
Наиболее важными среди условий однозначности являются краевые условия, разделяющиеся на пространственные (граничные) и временные (начальные). Их количество зависит от типа уравнений.
Дифференциальные уравнения в частных производных можно разделить на три типа:
I. Эллиптические уравнения описывают стационарные (установившиеся) процессы в замкнутой области пространства. Они не содержат производных по времени и к ним ставятся только граничные условия (на границах рассматриваемой области), которые в свою очередь могут быть разного рода. Наиболее важные:
а) Условия I рода, когда на границах области задаются значения искомых величин;
б) Условия II рода, когда на границах задаются значения производных от искомых величин;
в) Условия III рода, когда на границах области задаются соотношения между значениями искомых величин и их производных; например, соотношение
выражает так называемый «закон Ньютона-Рихмана» для теплообмена тела с омывающей его жидкостью температуры T0 (здесь a— постоянный коэффициент теплоотдачи).
Примером эллиптического уравнения может служить уравнение Лапласа.
II. Параболические уравнения описывают нестационарные процессы в замкнутой области. Они содержат первые производные искомых величин по времени и краевые условия к ним, помимо граничных, включают еще и начальные – значения искомых величин в начальный момент процесса. Примерами параболических уравнений служат уравнения теплопроводности и диффузии.
III. Гиперболические уравнения описывают волновые процессы. В них входят вторые производные по времени; соответственно, в состав краевых условий входят два начальных – значения искомых величин и их производных в начале процесса. Пример гиперболического уравнения – волновое уравнение.
Численные методы решения дифференциальных уравнений в частных производных.
Сложная структура дифференциальных уравнений в частных производных приводит к тому, что непосредственное решение их возможно только в редких случаях. Обычно прибегают к приближенным численным методам.
Рассмотрим в качестве простого примера стационарное распределение температуры в двумерной области (например, в тонкой длинной балке квадратного сечения, неравномерно подогреваемой и охлаждаемой с поверхности). Распределение температуры описывается дифференциальным уравнением теплопроводности, в котором отсутствуют производные по времени (процесс стационарный) и по координате z(в длинной балке можно не учитывать изменение температуры по длине). Поскольку l¹0, уравнение (44) приводится к виду
, (46)
т.е., к двумерному уравнению Лапласа.
Пусть справа и слева балка подогревается, а сверху и снизу охлаждается, например, по закону:
(47)
при x=±a,
(48)
приy=±a.
Выражения (47, 48) являются граничными условиями к уравнению (46).
Для численного решения дифференциальных уравнений в частных производных чаще всего пользуются различными вариантами метода сеток. Накроем рассматриваемую область сеткой из 2N+1горизонтальных и 2N+1вертикальных линий, равноотстоящих друг от друга. Расстояние между соседними линиями составит тогда h = a / N. Точки пересечения линий называются узлами сетки. Будем обозначать линии сетки номерами: i– для вертикальных (соответствующих определенному значению xi = ih, если номер iизменяется от –N до N) и j– горизонтальных ( соответствующих определенным значениям yj = jh). Тогда узел сетки можно определить по номерам линий, на пересечении которых он лежит.
Рассмотрим один из узлов сетки с номерами i,j.
Обозначим Ti,j=T(xi,yj). Аналогично обозначаются значения температуры или ее производных в других узлах сетки. Используя разложение функции T(x,y)в ряд Тейлора, можно записать:
Сложим эти выражения почленно:
.
Отсюда можно выразить вторую производную:
(49)
Такой способ получения выражений для производных сеточных функций является более общим, чем графический, представленный в п. 3.4.
Аналогично можно представить и вторую производную температуры по координатеy. Подставив приближенные выражения производных в уравнение (48), получим его разностное представление:
,
Дата добавления: 2016-12-16 ; просмотров: 1249 ; ЗАКАЗАТЬ НАПИСАНИЕ РАБОТЫ
Численные методы решения дифференциальных уравнений в частных производных
1.4.1. Использование метода сеток для решения параболических уравнений в частных производных
Одним из наиболее распространенных численных методов решения уравнений в частных производных является метод сеток. В методе сеток область Ω, в которой необходимо найти решение уравнения, разобьем прямыми, параллельными осям 
и 
, на прямоугольные области (рис.1.1), где
,
.
Точки, которые лежат на границе области Ω, называются внешними, остальные точки – внутренними. Совокупность всех точек называется сеткой, величины h и Δ – шагами сетки по хи t соответственно.
Рис.1.1. Сетка 
для области Ω
Идея метода сеток состоит в том, что вместо любой непрерывной функции 
будем рассматривать дискретную функцию 
, которая определена в узлах сетки , вместо производных функции будем рассматривать их простейшие разностные аппроксимации в узлах сетки. Таким образом, вместо системы дифференциальных уравнений в частных производных получим систему алгебраических уравнений. Чем меньше величины h иΔ, тем точнее получаемые алгебраические уравнения моделируют исходное дифференциальное уравнение в частных производных. В этом и последующих подразделах этой главы будет рассмотрен метод сеток для каждого из трех типов уравнений и его реализация в MATLAB. Знакомство с численными методами решения дифференциальных уравнений в частных производных начнем с разностных схем решения параболических уравнений.
Разностные схемы решения параболических уравнений будем рассматривать на примере следующего одномерного уравнения:
Построим сетку (см. рис.1.1). Для получения сеточного уравнения заменим производную 
приближенной разностной формулой:
.
В этой и последующих формулах 
– значение функции u в точке 
, 
– в точке 
, 
– в точке 
, 
– в точке 
, 
– в точке 
.
Для замены 
можно воспользоваться одной из приближенных разностных формул
;
.
Кроме того, заменим начальные и граничные условия их разностной аппроксимацией:
Заменив частные производные в задаче соответствующими разностными соотношениями, получим следующую вычислительную схему для расчета значений функции и в узлах сетки :
;
Это явная двухслойная разностная схема (рис.1.2).
Рис.1.2. Шаблон явной двухслойной разностной схемы
Учитывая, что на нулевом слое (при i = 0)все значения 
(как, впрочем, и 
) известны, можно сначала явно рассчитать значения 
, затем 
и так до 
.Для устойчивости разностной схемы значения шагов по t и х должны удовлетворять следующему условию:
.
Пример. Решить параболическое уравнение, описывающее распределение температуры в стержне длиной L,начальная температура стержня задается произвольной функцией 
. Температуры концов стержня равны 
.
Здесь а 2 – коэффициент температуропроводности; 
– коэффициент теплопроводности материала стержня; с – удельная теплоемкость; 
– плотность массы.
Подпрограмма решения задачи с помощью явной разностной схемы в MATLAB:
//Правая часть дифференциального уравнения
//Условие на левой границе
//Условие на правой границе
//Функция решения параболического уравнения методом сеток с помощью явной разностной схемы. N — количество участков разбиения интервала по х (0,1_); К – количество участков разбиения интервала по t (0,Т); а — коэффициент температуропроводности,
// Функция возвращает матрицу решений u вектора х, t
// Вычисляем шаг по х
// Вычисляем шаг по t
// Формируем массив х и первый столбец матрицы решений u
// из начального условия
//Формируем массив t, первую и последнюю строку матрицы решений u из граничных условий
gam = a ^ 2*delta/h ^ 2;
// Формируем матрицу решений u
Входными данными подпрограммы parabol решения задачи являются: N – количество интервалов, на которые разбивается отрезок (О, L); К – количество интервалов, на которые разбивается отрезок (О, Т); L – длина стержня, Т – интервал времени, а – параметр дифференциального уравнения. Функция возвращает три параметра: решение – сеточная функция u, определенная на сетке 
, массивы x и t.
График решения данной задачи изображен на рис.1.3.
Рис.1.3. График решения параболического уравнения при 
1.4.2. Использование метода сеток для решения гиперболических уравнений
Решение гиперболических уравнений также можно осуществить с помощью разностных схем. Разностные схемы решения одномерного гиперболического уравнения рассмотрим на примере следующего уравнения:
Построим сетку , в которой будем искать решение уравнения. Производную 
заменим соотношением:
,
остальные производные такие же, как и в предыдущем подразделе.
Получим следующую явную разностную схему решения уравнения:
Она устойчива при 
и по аналогии с разностной схемой (см. подразд. 1.4.1) может быть легко запрограммирована в MATLAB.
В качестве примера рассмотрим следующую задачу.
Пример. Решить начально-граничную задачу:
Ниже представлена функция giperbol решения данного уравнения, а на рис.1.4 – график полученного решения. Параметры функциианалогичны рассмотренным ранее, т.е. в подпрограмме решения параболического уравнения.
//Функция решения гиперболического уравнения с помощью явной разностной схемы. Входные данные: N — количество участков разбиения интервала (0,L); K — количество участков разбиения интервала (0,T); a — коэффициент теплопроводности; u — матрица решений в узлах сетки; массив x; массив t,
// Вычисляем шаг по х и по t
//Формируем массив x, первый и второй столбцы матрицы решений u из начального условия
//Формируем массив t, первую и последнюю строки матрицы решений u из граничных условий
//Формируем первую и последнюю строки матрицы решений u из граничных условий
//Формируем матрицу решений u
Рис.1.4. График решения гиперболического уравнения
1.4.3. Использование метода сеток для решения эллиптических уравнений
Разностные схемы решения одномерного эллиптического уравнения рассмотрим на следующем примере.
Пример. Рассмотрим разностную схему для эллиптического уравнения в прямоугольной области 
с граничными условиями Дирихле на границе Г:
Проведем в области Ω прямые, параллельные осям 
и 
, где 
; 
Для построения разностного уравнения заменим частные производные и граничные условия следующими соотношениями:
Преобразуем данную эллиптическую краевую задачу к следующей системе разностных уравнений:
Эту систему можно решать итерационными методами (например, методом Зейделя). В случае медленной сходимости итерационных процессов при решении сеточных уравнений, получаемых при аппроксимации гиперболических и эллиптических задач, можно попробовать заменить метод Зейделя градиентными методами (или методами релаксации). Ниже представлено решение эллиптического уравнения сеточным методом, а на рис.1.5 график полученного решения.
Рис.1.5. График решения эллиптического уравнения
// Функция ellip решения задачи
// Входные данные:R, a, b — значения, определяющие область решения задачи; Nx — количество участков разбиения интервала(R-b,R+b); Ny — количество участков разбиения интервала(-a,a); eps — точность решения уравнения методом Зейделя.
// psi — матрица решений в узлах сетки, массив x, массив y, k — количество итераций при решении разностного уравнения
// Вычисляем шаг по y и по x
// Формируем массив x, первый и последний столбцы матрицы решений psi из граничного условия
// Формируем массив y, первую и последнюю строки матрицы решений psi из граничного условия
// Вычисляем коэффициенты разностного уравнения
//Решение разностного уравнения методом Зейделя
http://lektsia.com/8x5c64.html