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

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

· для отрезка [a, b]:

(27)

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

3.5. Численное решение обыкновенных дифференциальных уравнений

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

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

При интегрировании ОДУ появляется ряд неопределенных констант интегрирования, число которых совпадает с порядком уравнения. Для определения этих констант служат начальные и граничные условия.

Начальные условия – это задание определенных связей для искомой функции при начальном значении аргумента (в начальный момент времени, в начальной точке какой-то линейно протяженной пространственной области). Задача отыскания решения системы ОДУ первого порядка с заданными начальными условиями называется в математике задачей Коши.

Граничные условия задаются при двух значениях аргумента. Обычно это условия на границах пространственной области.

а) тело падает с некоторой высоты. Необходимо определить зависимость скорости и положения тела над землей от времени.

На тело действуют две силы: сила тяжести, равная mg, где m – масса тела, g– ускорение свободного падения, и сила сопротивления воздуха, имеющая две составляющих – силу трения, пропорциональную скорости тела, а также вихревое сопротивление, пропорциональное квадрату скорости. Не вдаваясь в физические подробности, запишем общее сопротивление так: F_сопр=-(av+bv 2 ).Теперь уравнение второго закона Ньютона применительно к телу в проекции на ось Oy запишется (для удобства направим ось Оу вниз):

. (28)

. (29)

Таким образом, движение тела описывается уравнением (штрих, как обычно в математике, обозначает производную):

. (30)

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

. (31)

Задачу можно описать также системой двух уравнений (28), (29), каждое из которых имеет первый порядок. Таким образом, уравнение высокого порядка можно преобразовать в систему уравнений пониженного порядка, принимая производные за новые переменные аналогично (29). Начальные условия в этом случае будут выглядеть:

. (32)

Рассмотренная задача дает пример математической модели физического процесса, описываемого задачей Коши.

б) для улучшения охлаждения к корпусу теплового двигателя температурой Т₀ приварены ребра – тонкие металлические пластины. Необходимо определить распределение температуры вдоль ребра и количество тепла, отводимое с его помощью в окружающую среду.

Распространение тепла в толще ребра описывается так называемым законом Фурье:

, (33)

где Q – количество тепла, переносимое за 1 сек через сечение ребра площадью F=bH; l —так называемый коэффициент теплопроводности, Т – температура. Количество тепла, отводимое в окружающую среду с участка поверхности ребра длиной dx, равно (Т_ос – температура окружающей среды):

.

Первое слагаемое в квадратных скобках выражает охлаждение поверхности ребра омывающим его воздухом, второе – отвод тепла путем теплового излучения. Исходя из теплового баланса слоя ребра, ограниченного сечениями, соответствующими значениям координат x и x+dx, можем записать:

,

В левой части равенства стоит выражение производной

(так как наличие дифференциала dxуказывает на предельный переход при х®0, ср. (22)). С учетом (33) получаем уравнение:

. (34)

Это уравнение второго порядка и для него надо поставить два условия. Одно из них очевидно: T=T₀при х=0. Второе условие в теории теплообмена ставится по-разному. Простейший вариант (пригодный для длинных ребер) – считать, что температура конца ребра совпадает с температурой окружающей среды: Т=Т_оспри х=L. Таким образом, в данной задаче для дифференциального уравнения (34) поставлены граничные условия.

Рассмотрим некоторые методы решения ОДУ. Начнем с методов решения задачи Коши, то есть, численного решения дифференциальных уравнений вида

(35)

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

(36)

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

или (37)

Используя формулу производной сеточной функции (23), можно записать:

Это и есть расчетная формула метода Эйлера. Для решения отрезок оси Ох,на котором надо найти решение уравнения, разбивают на малые шаги h (создают сетку), затем в правую часть подставляют значения и . Это позволяет приближенно определить значение искомой функции y(x) в точке . Затем найденные значения и опять подставляют в правую часть формулы (6), вычисляют и т.д. Если обозначить , то формулу (6) можно записать в окончательном виде

(38)

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

(39)

(40)

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

Метод Эйлера реализуется существенно легче, чем метод Рунге-Кутта, поэтому его можно применить даже при расчете вручную, однако его точность крайне низка – он имеет первый порядок точности, тогда как метод Рунге-Кутта – четвертый. Поэтому в компьютерных расчетах используется метод Рунге-Кутта или другие методы высокого порядка точности.

Для решения краевых задач применяются разные методы. Рассмотрим самый простой из них – метод пристрелки или стрельбы. Он основан на преобразовании краевой задачи в задачу Коши.

Рассмотрим для примера задачу (34), которую запишем в виде:

(41)

Вместо последнего граничного условия при х=L зададим дополнительное условие при х=0: q = q₀. Тогда мы получим задачу Коши. Величину q₀ можно выбрать произвольно, например, принять в первом приближении равной нулю.

Решив задачу Коши одним из вышеописанных методов, мы определим значение неизвестной величины Т на конце ребра, при х=L. Скорее всего, оно будет отличаться от значения Т_ос. Сущность метода пристрелки заключается в том, что значение дополнительного начального условия подбирается таким образом, чтобы после решения задачи Коши оказалось выполненным граничное условие на конце рассматриваемого отрезка (в данном случае при x=L). Это действительно напоминает пристрелку при артиллерийской стрельбе, где дополнительному начальному условию можно сопоставить угол наклона ствола орудия, который меняется до тех пор, пока не будет выполнено граничное условие на другом конце – совпадение траектории полета снаряда с точкой местонахождения цели.

Фактически получается, что значение граничного условия на конце отрезка является функцией величины дополнительного начального условия: T½_x=L = f(q₀) . Тогда условие совпадения конечного условия с заданной величиной Т_ос запишется в виде некоторого нелинейного уравнения:

,

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

3.6. Моделирование процессов, приводящих
к дифференциальным уравнениям в частных производных

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

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

(42)

Такое непрерывное распределение температуры по точкам в сплошной среде называется полем температуры. Поля могут быть образованы и другими величинами – давлением, концентрацией растворённого вещества и т.п.

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

.

Такое векторное поле можно представить, как комбинацию трех скалярных полей проекций скорости на оси координат.

3.7. Дифференциальное уравнение теплопроводности

Рассмотрим процесс теплообмена в твердом теле. Количество тепла, передаваемое за некоторое бесконечно малое время dtчерез бесконечно малый элемент dF поверхности, нормальной к координатной оси Ox, пропорционально dtи dF, а также зависит от скорости изменения температуры в направлении оси Ox– чем быстрее изменяется температура, тем больше передается тепла. В большинстве случаев связь количества передаваемого тепла со скоростью изменения температуры вдоль Ox принимается линейной:

. (43)

Коэффициент пропорциональности l называется коэффициентом теплопроводности; знак «минус» указывает на известный из физики факт – тепло передается в сторону убывания температуры.

называется плотностью потока тепла в направлении координаты x (численно величина q_xсоответствует количеству тепла, передаваемому за 1 сек через 1 м 2 поверхности). Аналогично можно записать:

Величины q_x, q_y, q_z можно рассматривать, как компоненты вектора

,

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

.

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

В векторной форме закон теплопроводности записывается:

.

Это соотношение называется законом Фурье. Аналогичный вид имеет так называемый закон Фика, описывающий диффузию какого-либо вещества в сплошной среде:

,

где вектор— вектор плотности потока массы i—компонента смеси веществ (модуль этого вектора показывает, сколько данного вещества проходит за 1 сек через малый элемент поверхности, нормальной к направлению переноса, в расчете на 1 м 2 площади поверхности), grad c_i– градиент концентрации данного вещества, r— плотность среды, D_i – коэффициент диффузии данного вещества.

В силу аналогии законов Фурье и Фика многие закономерности и расчетные формулы теплообмена оказываются справедливыми и для массообмена (диффузии).

Рассмотрим бесконечно малый элемент сплошной среды в виде параллелепипеда с ребрами dx, dy, dz. Через левую грань объема в него поступает количество теплоты dQ_x, описываемое формулой (43). Одновременно через правую грань из него уходит теплота

Если температура меняется вдоль оси Ох, значения производных на противоположных гранях неодинаковы. В итоге внутри рассматриваемого элемента остается тепло

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

.

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

.

Накопление тепла вызывает нагрев вещества в объёме:

,

где — частный дифференциал температуры.

Приравняв два выражения для d 2 Q, получаем

(44)

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

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

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

;

— волновое уравнение (уравнение распространения волн в сплошной среде):

,

где u – величина, волновое колебание которой рассматривается (давление, температура, электрический потенциал и др.), v– скорость распространения волны;

, (45)

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

3.8. Краевые задачи для уравнений в частных производных

Уравнение теплопроводности применимо к анализу распространения тепла в любой системе – от куска сыра до стенки атомного реактора. Поэтому при рассмотрении конкретных задач к уравнению надо добавить описание особенностей задачи, отличающих ее постановку от других задач теплопроводности.

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

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

Дифференциальные уравнения в частных производных можно разделить на три типа:

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

а) Условия I рода, когда на границах области задаются значения искомых величин;

б) Условия II рода, когда на границах задаются значения производных от искомых величин;

в) Условия III рода, когда на границах области задаются соотношения между значениями искомых величин и их производных; например, соотношение

выражает так называемый «закон Ньютона-Рихмана» для теплообмена тела с омывающей его жидкостью температуры T₀ (здесь 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– для вертикальных (соответствующих определенному значению x_i = ih, если номер iизменяется от –N до N) и j– горизонтальных ( соответствующих определенным значениям y_j = jh). Тогда узел сетки можно определить по номерам линий, на пересечении которых он лежит.

Рассмотрим один из узлов сетки с номерами i,j.

Обозначим T_i,j=T(x_i,y_j). Аналогично обозначаются значения температуры или ее производных в других узлах сетки. Используя разложение функции T(x,y)в ряд Тейлора, можно записать:

Сложим эти выражения почленно:

.

Отсюда можно выразить вторую производную:

(49)

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

Аналогично можно представить и вторую производную температуры по координатеy. Подставив приближенные выражения производных в уравнение (48), получим его разностное представление:

,

Дата добавления: 2016-12-16 ; просмотров: 1248 ; ЗАКАЗАТЬ НАПИСАНИЕ РАБОТЫ

Численное решение уравнений в частных производных эллиптического типа на примере уравнений Лапласа и Пуассона

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

6.1. Постановка задачи. Простейшая разностная схема «крест». Устойчивость схемы «крест»

Будем рассматривать двухмерное уравнение Пуассона

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

( — заданная на границе функция ).

В случае прямоугольной области граничные условия удобно записать в следующем виде:

Для простоты выкладок введем равномерную расчетную сетку с узлами m, y_l> , m, l = 0, 1, . , M с равным количеством шагов по каждому пространственному направлению, сеточную область D — совокупность всех узлов сетки, включая граничные, и сеточную функцию < u_ml >. В этом случае шаги по координатам предполагаются равными. В случае неравных шагов по каждому направлению полученные результаты не изменятся, а запись уравнений станет более громоздкой.

Выбираем простейший пятиточечный шаблон разностной схемы «крест» . На этом шаблоне аппроксимирующее разностное уравнение легко выписать. Для этого производные заменим вторыми разностями:

где h — шаг по координатам, или в операторной форме

Эту же разностную схему можно записать в каноническом виде для разностных схем для эллиптических уравнений:

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

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

Здесь учтено разложение проекции точного решения в ряд Тейлора

и аналогичное разложение для u_{m — 1}.

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

Рассмотрим устойчивость полученной схемы. Отметим, что методы исследования на устойчивость , применяемые для эволюционных (зависящих от времени) уравнений, здесь не работают. Действовать приходится на основе определения устойчивости.

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