Приведем в разностной схеме (7) подобные слагаемые, перенеся в правую часть значения сеточной функции с индексом k (как часто говорят, с предыдущего слоя по времени), а в левую — с индексом k+i (т. е. со следующего временного слоя). Кроме этого, введем коэффициент с, который будет характеризовать отношение шагов разностной схемы по времени и пространству
Несколько забегая вперед, заметим, что значение параметра с, называемого коэффициентом Куранта, имеет большое значение для анализа устойчивости разностной схемы. С учетом этих замечаний, разностная схема (7) запишется в виде:
(8)Рис. 6. Шаблон аппроксимации явной схемы для уравнения теплопроводности
Множители для каждого из значений сеточной функции в узлах шаблона, соответствующие разностному уравнению (8), приведены рядом с каждой точкой шаблона на рис. 6. Фактически геометрия шаблона и эти множители задают построенную нами разностную схему.
Несложно убедиться в том, что для получения замкнутой системы разностных алгебраических уравнений систему (8) необходимо дополнить дискретным представлением начального и граничных условий. Тогда число неизвестных будет в точности равно числу уравнений, и процесс формирования разностной схемы будет окончательно завершен.
Если присмотреться к разностным уравнениям (8) повнимательнее, то можно сразу предложить несложный алгоритм реализации этой разностной схемы. Действительно, каждое неизвестное значение сеточной функции со следующего временного слоя, т. е. левая часть соотношения (8) явно выражается через три ее значения с предыдущего слоя (правая часть), которые уже известны. Таким образом, в случае уравнения теплопроводности нам очень повезло. Для расчета первого слоя по времени следует попросту подставить в (8) начальное условие (известные значения и с нулевого слоя в узлах сетки), для расчета второго слоя достаточно использовать вычисленный таким образом набор и с первого слоя и т. д. Из-за того что разностная схема сводится к такой явной подстановке, ее и называют явной, а благодаря пересчету значений с текущего слоя через ранее вычисленные слои — схемой бегущего счета.
Линейное уравнение
Сделанные замечания относительно реализации явной схемы для уравнения диффузии тепла сразу определяют алгоритм ее программирования в Mathcad. Для решения задачи нужно аккуратно ввести в листинг соответствующие формулы при помощи элементов программирования.
Решение системы разностных уравнений (8) для модели без источников тепла, т.е. ф(x,T,t)=0 и постоянного коэффициента диффузии D=const приведено в листинге 1. В его первых трех строках заданы шаги по временной и пространственной переменным t и А, а также коэффициент диффузии о, равный единице. В следующих двух строках заданы начальные (нагретый центр области) и граничные (постоянная температура на краях) условия соответственно. Затем приводится возможное программное решение разностной схемы, причем функция пользователя v(t) задает вектор распределения искомой температуры в каждый момент времени (иными словами, на каждом слое), задаваемый целым числом t.
Начальное распределение температуры вдоль расчетной области и решение для двух моментов времени показано на рис. 7 сплошной, пунктирной и штриховой линиями соответственно. Физически такое поведение вполне естественно — с течением времени тепло из более нагретой области перетекает в менее нагретую, а зона изначально высокой температуры остывает и размывается.
Листинг 11. Явная схема для линейного уравнения теплопроводности
Рис. 7. Решение линейного уравнения теплопроводности (продолжение листинга 1)
Нелинейное уравнение
Намного более интересные решения можно получить для нелинейного уравнения теплопроводности, например, с нелинейным источником тепла ф(u)=103(u-u3). Заметим, что в листинге 1 мы предусмотрительно определили коэффициент диффузии и источник тепла в виде пользовательских функций, зависящих от аргумента и, т. е. от температуры. Если бы мы собирались моделировать явную зависимость их от координат, то следовало бы ввести в пользовательскую функцию в качестве аргумента переменную х, как это сделано для источника тепла ф. Поэтому нет ничего проще замены определения этих функций с констант D(U)=1 и ф(х,u)=о на новые функции, которые станут описывать другие модели диффузии тепла. Начнем с того, что поменяем четвертую строку листинга 1 на ф(х,и)=ю3-(и-и3), не изменяя пока постоянного значения коэффициента диффузии.
Рис. 7. Решение уравнения теплопроводности с нелинейным источником (тепловой фронт)
Рис. 8. Решение уравнения теплопроводности с нелинейным источником и коэффициентом диффузии (режим локализации горения)
Читателю предлагается поэкспериментировать с этим и другими нелинейными вариантами уравнения теплопроводности. Существенно, что такие интересные результаты удается получить лишь численно, а в Mathcad только с применением элементов программирования.
В отличие от явной схемы Эйлера, неявная является безусловно-устойчивой (т.е. не выдающей "разболтки" ни при каких значениях коэффициента Куранта). Однако ценой устойчивости является необходимость решения на каждом шаге по времени системы алгебраических уравнений.
Построение неявной разностной схемы
Чтобы построить неявную разностную схему для уравнения диффузии, используем шаблон, изображенный на рис. 9, т. е. для дискретизации пространственной производной будем брать значения сеточной функции с верхнего (неизвестного) слоя по времени. Таким образом, разностное уравнение для (i,k)-ro узла будет отличаться от уравнения для явной схемы (7) только индексами по временной координате в правой части:
(9)Если привести подобные слагаемые, то получится система уравнений, связывающая для каждого 1-го узла три неизвестных значения сеточной функции (в самом этом узле и в соседних с ним слева и справа узлах). Множители при неизвестных значениях сеточной функции в узах шаблона показаны на рис. 9 в виде подписей, подобно тому, как это было сделано для явной схемы.
Рис. 9. Шаблон неявной схемы для уравнения теплопроводности
Очень важно, что если само уравнение теплопроводности линейно, то с в левой части разностного уравнения является константой, а ф в его правой части может зависеть только от первой степени и. Поэтому система уравнений (10) для всех пространственных узлов 1=1. .м-l является линейной системой, что существенно упрощает ее решение (поскольку известно, что для линейных систем с ненулевым определителем решение существует и является единственным). Напомним, что для получения замкнутой системы линейных уравнений необходимо дополнить данный набор разностных уравнений граничными условиями, т.е. известными значениями сеточной функции для i=0 и i=M.
Рис. 9. Решение линейного уравнения теплопроводности при помощи неявной схемы на первом слое по времени (листинг 2)
Листинг 2. Неявная схема для линейного уравнения теплопроводности
Для реализации неявной схемы, таким образом, можно использовать комбинацию средств программирования Mathcad и встроенной функции решения системы линейных уравнений isolve. Один из возможных способов решения предложен в листинге 2. Большая часть этого листинга является вводом параметров задачи (шагов, начальных и граничных условий), и только в последней его строке определяется функция пользователя, вычисляющая сеточную функцию на каждом временном слое (при помощи встроенной функции решения системы линейных уравнений isolve). В нескольких предыдущих строках листинга (после расчета коэффициента Куранта) формируется матрица системы уравнений, которая записывается в подходящем для Mathcad виде, как это сделано в листинге 2. Как несложно убедиться, столбец правых частей разностных уравнений выражается вычисленными значениями сеточной функции с предыдущего слоя.
В данной работе приведены некоторые примеры применения дифференциальных уравнений для моделирования таких реальных процессов, как колебания струны, электрические колебания в проводах, распространение тепла в стержне и пространстве, распространение температурных волн в почве, дифракция излучения на сферической частице.
Во многих областях физики, математики и других естественных наук часто используются численные и эмпирические методы для решения прямых и обратных задач. Следует отметить особую роль дифференциальных уравнений при решении таких задач, поскольку не всегда удается установить функциональную зависимость между искомыми и данными переменными величинами, но зато часто удается вывести дифференциальное уравнение, позволяющее точно предсказать протекание определенного процесса при определенных условиях.