Кирьянов. Самоучитель MathCad 11. Кирьянов д в
 Скачать 10.75 Mb.
|
|
{ х- -1 1 i ± + 1 Шаблон неявной схемы для уравнения теплопроводности 332 Часть III. Численные методы Если привести подобные слагаемые, то получится система уравнений, связывающая для каждого узла три неизвестных значения сеточной функции (в самом этом узле ив соседних с ним слева и справа узлах. Множители при неизвестных значениях сеточной функции в узлах шаблона показаны на рис. 13.11 в виде подписей, подобно тому как это было сделано для явной схемы (см. Очень важно, что если само уравнение теплопроводности линейно, то св левой части разностного уравнения является константой, а ф в его правой части может зависеть только от первой степени и. Поэтому система уравнений) для всех пространственных узлов является линейной системой, что существенно упрощает ее решение (поскольку известно, что для линейных систем с ненулевым определителем решение существует и является единственным. Напомним, что для получения замкнутой системы линейных уравнений необходимо дополнить данный набор разностных уравнений граничными условиями, те. известными значениями сеточной функции и Примечание Если рассматривать нелинейный случай, тона каждом шаге повремени пришлось бы решать систему нелинейных уравнений, число решений которых могло бы быть большими среди них требовалось бы отыскать нужное, а не паразитное решение. Для реализации неявной схемы, таким образом, можно использовать комбинацию средств программирования и встроенной функции решения системы линейных уравнений Один из возможных способов решения предложен в листинге 13.2. Большая часть этого листинга является вводом параметров задачи (шагов, начальных и граничных условий, и только в последней его строке определяется функция пользователя, вычисляющая сеточную функцию на каждом временном слое (при помощи встроенной функции решения системы линейных уравнений lsoive). В нескольких предыдущих строках листинга (после расчета коэффициента Куранта сои) формируется матрица системы уравнений А, которая записывается в подходящем для Mathcad виде, как это сделано в листинге. Как несложно убедиться, столбец правых частей разностных уравнений выражается вычисленными значениями сеточной функции с предыдущего слоя. Результаты расчетов по неявной схеме показаны на рис. 13.12 и, как видно, они дают примерно те же результаты, что ив случае применения явной схемы (см. рис. 13.7). Обратите внимание, что решение устойчиво при любых значениях коэффициента Куранта (в том числе, и больших 1), поскольку, как следует из соответствующих положений теории численных методов, неявная схема является безусловно устойчивой Глава 13. Уравнения в частных производных u) :=0 D:=l Border :=0 (х) т:= 0.005 ДМ Сои , 0 V ( n Сои u i f n = О - Алгоритм прогонки Приведем в данном разделе описание чрезвычайно популярного алгоритма реализации неявных разностных схем, который называется методом прогонки Этот алгоритм имел историческое значение для становления технологий расчетов уравнений в частных производных, и мы просто не можем не упомянуть о нем в этой книге Примечание Сразу оговоримся, что его применение для решения уравнений в частных производных в среде Mathcad может быть оправдано, только если Выработаете сочень частыми сетками, которые приводят к системам разностных уравнений большой размерности и, соответственно, очень долгому времени вычислений. Основным вычислительным ядром программы, реализующей на Mathcad неявную разностную схему, было решение (на каждом временном слое) системы линейных алгебраических уравнений, задаваемых матрицей АЗа- метим, что эта матрица, как говорят, имеет диагональное преобладание а Часть III. Численные методы точнее, является (рис. 13.13). Все ее элементы, кроме элементов на главной диагонали и двух соседних диагоналях, равны нулю. Сточки зрения оптимизации быстродействия алгоритма, применение встроенной функции является весьма расточительным, поскольку основной объем арифметических операций, выполняемых компьютером (а он составляет, как нетрудно убедиться величину порядка сводится к непроизводительному перемножению нулей i г 13.12. Решение линейного уравнения теплопроводности при помощи неявной схемы на первом слое повремени (листинг Для отыскания решения линейных систем алгебраических уравнений имеется чрезвычайно эффективный алгоритм, называемый прогонкой, который позволяет снизить число арифметических операций на целый порядок, т. е. до значения порядка Это означает, что при использовании пространственных сеток с юоо узлами выигрыш во времени вычислений составит величину порядка ю Реализация данного алгоритма приведена в листинге, который является продолжением листинга 13.2, используя определенные в нем коэффициенты матрицы А, а также начальное условие. Программа листинга 13.3 осуществляет пересчет одного шага повремените. заменяет содержимое столбца с предыдущего временного слоя вычисленными значениями неизвестной функции со следующего слоя. Первые пять строк листинга 13.3 представляют так называемый обратный ход прогонки, а последние две строки — ее прямой ход. Заинтересовавшемуся читателю предлагается самому оформить представленный алгоритм прогонки в виде программы решения разностных уравнений для вычисления произвольного временного слоя по примеру листингов 13.1 и 13.2. Заметим, что описание этого знаменитого алгоритма можно отыскать практически в любом современном учебнике по численным методам Глава 13. Уравнения в частных производных 0 j 0 •? 11 fl -2 0 5 о 0 0 о Рис. 13.13. Матрица системы линейных разностных уравнений для неявной схемы (листинг 13.2 для Листинг 13.3. Алгоритм прогонки (продолжение листинга := т т : = 0 . . 13.2.3. О возможности решения многомерных уравнений Все, что было сказано до сих пор, касалось исключительно способов решения одномерных (в смысле пространственных координат) уравнений. И алгоритмы разностных схем, и встроенные функции, включая появившиеся в 11-й версии (см. следующий разд относились к уравнениям, зависящим от одной пространственной координаты. Можно ли при помощи Mathcad решать двумерные или трехмерные (пространственные) уравнения Сточки зрения программирования пользователем численных алгоритмов типа метода сеток, принципиальных ограничений нет. Разумеется, если сначала аккуратно выписать разностную схему 336 Часть Численные методы соответствующего многомерного дифференциального уравнения, то вполне возможно запрограммировать ее при помощи описанных нами средств. Самым главным противодействием будет существенное увеличение времени расчетов. Простая оценка необходимого количества операций показывает, что ввод зависимости уравнения от второй пространственной координаты многократно увеличивает число разностных уравнений, которые должны решаться при реализации каждого шага повремени. К примеру, если используется пространственная сетка из 100 узлов по каждой координате, то вместо 10 2 разностных уравнений на каждом шаге придется решать уже 10 уравнений, те. объем вычислений сразу же возрастает враз. Вообще говоря, пакет Mathcad не является экономичной средой вычислений, ибо- роться сих сильно возрастающим объемом пользователю следует еще на этапе разработки алгоритма. Хорошим примером такой борьбы может служить применение специфических алгоритмов, типа метода прогонки (см. разд. Алгоритм прогонки этой главы). Приведем некоторые дополнительные замечания, связанные с возможностью осуществить редукцию громоздких (в смысле организации вычислений) двумерных задач к более простым. Рассмотрим ради примера двумерное уравнение теплопроводности (1) без источника с нулевыми граничными условиями и некоторым начальным двумерным распределением температуры по расчетной поверхности у, t) _ у, t) у, Произведем дискретизацию данного уравнения повременной координате, заменяя первую производную ее разностным аналогом и несколько перегруппировывая слагаемые и множители - Эх ду 2 Как Вы видите, мы используем неявную разностную схему, заранее заботясь о том, чтобы разностная задача была более устойчивой. Здесь — известная с предыдущего повремени функция двух пространственных переменных, a функция, подлежащая определению при реализации каждого шага по времени. Можно посмотреть на полученную задачу с несколько другой стороны — а именно как на дифференциальное уравнение относительно неизвестной функции двух переменных Подчеркнем, что такое уравнение получается для каждого шага повремените. для реализации всей разностной схемы требуется решить большое число таких уравнений. С предложенной точки зрения, на каждом временном шаге необходимо решить некоторое двумерное эллиптическое линейное уравнение, причем его граничные условия определяются граничными условиями исходной задачи Глава 13. Уравнения в частных производных Это уравнение очень похоже на уравнение Пуассона стой лишь разницей, что в его правую часть, описывающую источник, входит неизвестная функция (к счастью, линейно. Таким образом, зависимость от найденного на предыдущем шаге повремени решения определяется зависимостью от него источника, те. правой части. Суммируя сказанное, можно констатировать, что если Вы имеете запрограммированный алгоритм решения выписанного эллиптического уравнения, чуть более сложного, чем уравнение Пуассона, то его с легкостью можно использовать в качестве подпрограммы реализации разностной схемы двумерного уравнения теплопроводности. Забегая вперед, приходится отметить, что, к сожалению, встроенные функции Mathcad для решения уравнения Пуассона в данном случае не годятся в качестве такой подпрограммы, поскольку предполагают независимость источника от самой неизвестной функции и могут справляться лишь с правой частью, зависящей только от пространственных координат (см. разд. Не будем далее останавливаться на способах решения многомерных уравнений, ограничившись этими замечаниями относительно путей оптимизации алгоритмов их решения Встроенные функции для решения уравнений в частных производных Как видно из предыдущего раздела, с уравнениями в частных производных вполне можно справиться, и не прибегая к специфическим средствам. Между тем, в пакете Mathcad 11 имеется несколько встроенных функций, при помощи которых можно автоматизировать процесс решения дифференциальных уравнений в частных производных. Рассмотрим в данном разделе основные аспекты их применения, отмечая не только инструкции по их применению, описанные разработчиками Mathcad, но также и некоторые скрытые возможности этих функций. Параболические и гиперболические уравнения В новой версии Mathcad 11 разработчики впервые применили встроенную функцию для решения уравнений в частных производных, отлично осознавая значимость этих задач для современного исследователя и инженера. Эта функция применяется в рамках вычислительного блока, начинающегося ключевым словом Given и пригодна для решения различных гиперболических и параболических уравнений. Встроенная функция для решения одномерного уравнения (или системы уравнений) в частных производных (того, которое определит пользователь в Часть ill. Численные методы рамках вычислительного блока Given), зависящего от времени t и пространственной координаты х, имеет целый набор различных аргументов и работает следующим образом Pdesolve(u, x, xrange, t, trange, [xpts] , [tpts] ) ) — возвращает скалярную (для единственного исходного уравнения) или векторную (для системы уравнений) функцию двух аргументов (x,t), являющуюся решением дифференциального уравнения (или системы уравнений) в частных производных. Результирующая функция получается интерполяцией сеточной функции, вычисляемой согласно разностной схеме. • и — явно заданный вектор имен функций (без указания имен аргументов, подлежащих вычислению. Эти функции, а также граничные условия (в форме Дирихле или Неймана) должны быть определены пользователем перед применением функции в вычислительном блоке после ключевого слова Given. Если решается не система уравнений в частных производных, а единственное уравнение, то, соответственно, вектор должен содержать только одно имя функции и вырождается в скаляр. • х пространственная координата (имя аргумента неизвестной функции пространственный интервал, те. вектор значений аргументах для граничных условий. Этот вектор должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала — время (имя аргумента неизвестной функции расчетная временная область вектор значений аргумента который должен состоять из двух действительных чисел (представляющих левую и правую границу расчетного интервала повремени количество пространственных точек дискретизации (может не указываться явно, в таком случае будет подобрано программой автоматически количество временных слоев, те. интервалов дискретизации повремени (также может не указываться пользователем явно). Примечание Помимо этой функции для решения параболических и гиперболических уравнений, начиная с новой версии 11, можно использовать еще одну встроенную функцию ( ) . Функция () имеет еще большее число аргументов и позволяет управлять дополнительными параметрами метода сеток. Однако пользоваться ею намного сложнее, чем функцией P d e s o l v e O , и поэтому в нашей книге мы не будем на ней особо останавливаться Глава 13. Уравнения в частных производных В качестве примера использования этой новой функции Mathcad 11 листинг) используем тоже самое одномерное уравнение теплопроводности) с граничными и начальными условиями (6) и (7). | Листинг 13.4. Решение одномерного уравнения теплопроводности | D : = 0 . 1 G i v e n ( x , t ) • t ) u ( x , 0 ) Ф { x - 0 . 4 5 ) - ф ( х - 0 . 5 5 ) u ( 0 , t ) • 0 u ( L , t ) = 0 u | . , Для корректного использования функции предварительно, после ключевого слова Given, следует записать само уравнение и граничные условия при помощи логических операторов (для их ввода в Mathcad существует специальная панель. Обратите внимание, что уравнение должно содержать имя неизвестной функции вместе с именами аргументов (а не так, как она записывается в пределах встроенной функции Для идентификации частных производных в пределах вычислительного блока следует использовать нижние индексы, например для обозначения второй производной функции по пространственной координате х.. Как видно из рис. 13.14, на котором изображены результаты расчетов по листингу 13.4, встроенная функция с успехом справляется с уравнением диффузии, отыскивая уже хорошо знакомое нам решение. Заметим, что использование встроенной функции pdesolve связано с довольно громоздкими вычислениями, которые могут отнимать существенное время. Примечание Как Вы можете заметить, выбирать величину шага по пространственной и временной переменным может как сам алгоритм, таки пользователь (неявным образом, через число узлов сетки. Читателю предлагается повторить вычисления листинга 13.4 для различных комбинаций параметров (главным образом, числа узлов сетки, чтобы проверить, в каких случаях алгоритм встроенной функции справляется с задачей, выдавая верное решение, а в каких дает сбой. Приведем еще один пример применения функции Pdesolve для решения уравнений в частных производных. Рассмотрим одномерное волновое уравнение, которое описывает, например, свободные колебания струны музыкального инструмента t) t) . (11) Часть III. Численные методы Рис. 13.14. Решение уравнения диффузии тепла при помощи встроенной функции Pdesolve (листинг Здесь неизвестная функция u(x,t) описывает динамику смещения профиля струны относительно невозмущенного (прямолинейного) положения, а параметр с характеризует материал, из которого изготовлена струна. Как Вы видите, уравнение (И) содержит производные второго порядка как по пространственной координате, таки повремени. Для того чтобы можно было использовать встроенную функцию pdesolve, необходимо переписать волновое уравнение в виде системы двух уравнений в частных производных, введя вторую неизвестную функцию Программа для решения волнового уравнения приведена в листинге 13.5, а результат — на рис. Листинг 13.5. Решение волнового уравнения с:=1 : = 2 л ( x , t ) с хи Глава Уравнения в частных производных Т Рис. 13.15. Решение волнового уравнения (листинг 13.5) 13.3.2. Эллиптические уравнения Решение эллиптических уравнений в частных производных реализовано только для единственного типа задач — двумерного уравнения Пуассона. Это уравнение содержит вторые производные функции по двум пространственным переменным: (12) Эх , — + у) Уравнение Пуассона описывает, например, распределение электростатического поля и(х,у) в двумерной области с плотностью заряда f(x,y) или (см. разд стационарное распределение температуры и(х,у) на плоскости, в которой имеются источники (или поглотители) тепла с интенсивностью Примечание Несмотря на то, что применение встроенных функций, описанных в данном разделе, анонсировано разработчиками Mathcad только для уравнения Пуассона, их можно применять и для решения других уравнений, даже необязательно эллиптического типа. О том, как осуществить такие расчеты, написано в конце данного раздела 342 Часть III. Численные методы Уравнение Пуассона с нулевыми граничными условиями Корректная постановка краевой задачи для уравнения Пуассона требует задания граничных условий. В решение ищется на плоской квадратной области, состоящей из точек. Поэтому граничные условия должны быть определены пользователем для всех четырех сторон упомянутого квадрата. Самый простой вариант — это нулевые граничные условия, т. е. постоянная температура по всему периметру расчетной области. В таком случае можно использовать встроенную функцию — матрица решения уравнения Пуассона размера на квадратной области с нулевыми граничными условиями матрица размера задающая правую часть уравнения Пуассона; • — параметр численного алгоритма (количество циклов в пределах каждой итерации). Внимание! Сторона квадрата расчетной области должна включать точно шагов, те узлов, где п — целое число. Параметр численного метода ncycie в большинстве случаев достаточно взять равным 2. Листинг 13.6 содержит пример использования функции muitigrid для расчета краевой задачи на области ззхзз точки и точечным источником тепла вместе, задаваемом координатами (15,20) внутри этой области. |