Кирьянов Д. - MathCad 11 (1077323), страница 54
Текст из файла (страница 54)
13.10. Численное решение уравнения теплопроводности при помощи явнойсхемы Эйлера (см. листинг 13.1 ниже с временным шагом т=0.0015)Глава 13. Уравнения в частных производных331Характерная «разболтка» решения как раз и является проявлением неустойчивости явной схемы Эйлера для выбранного соотношения шагов по времени и пространству. В теории численных методов показывается, что явнаясхема Эйлера для уравнения теплопроводности устойчива при значенияхкоэффициента Куранта, меньших 1, и неустойчива в противоположном случае.
Иными словами, существует ограничение для выбора соотношения шагов, заключающееся в том, что для расчета на более частых пространственных сетках необходимо использовать также и малые шаги по времени.ПримечаниеКак несложно убедиться, для т=0.0005 коэффициент Куранта с = 0 . 4 , длят=0.0010 он все еще меньше единицы: с=0. 8, а для т=0.0015 решение ужебольше единицы: с=1.2, в связи с чем схема становится неустойчивой(рис. 13.10).13.2.2. Неявная схема ЭйлераВ отличие от явной схемы Эйлера, неявная является безусловно устойчивой(т. е.
не выдающей «разболтки» ни при каких значениях коэффициента Куранта). Однако, ценой устойчивости является необходимость решения накаждом шаге по времени системы алгебраических уравнений.Построение неявной разностной схемыЧтобы построить неявную разностную схему для уравнения диффузии, используем шаблон, изображенный на рис.
13.11, т. е. для дискретизации пространственной производной будем брать значения сеточной функции сверхнего (неизвестного) слоя по времени. Таким образом, разностное уравнение для (i,k)-ro узла будет отличаться от уравнения для явной схемы (8)только индексами по временной координате в правой части:" » « - " » .„;.-»« - 2 uwк+1u,(W1+{,(101-U1*кх--1i±+1Рис. 13.11. Шаблон неявной схемыдля уравнения теплопроводности)332Часть III. Численные методыЕсли привести подобные слагаемые, то получится система уравнений, связывающая для каждого i-ro узла три неизвестных значения сеточной функции (в самом этом узле и в соседних с ним слева и справа узлах). Множители при неизвестных значениях сеточной функции в узлах шаблонапоказаны на рис. 13.11 в виде подписей, подобно тому как это было сделано для явной схемы (см. рис.13.6).Очень важно, что если само уравнение теплопроводности линейно, то с влевой части разностного уравнения является константой, а ф в его правойчасти может зависеть только от первой степени и.
Поэтому система уравнений (10) для всех пространственных узлов i = i . . . M - i является линейнойсистемой, что существенно упрощает ее решение (поскольку известно, чтодля линейных систем с ненулевым определителем решение существует иявляется единственным). Напомним, что для получения замкнутой системылинейных уравнений необходимо дополнить данный набор разностныхуравнений граничными условиями, т.
е. известными значениями сеточнойфункции для i=o и i=M.ПримечаниеЕсли рассматривать нелинейный случай, то на каждом шаге по времени пришлось бы решать систему нелинейных уравнений, число решений которых могло бы быть большим, и среди них требовалось бы отыскать нужное, а не паразитное решение.Для реализации неявной схемы, таким образом, можно использовать комбинацию средств программирования Mathcad и встроенной функции решения системы линейных уравнений lsoive.
Один из возможных способов решения предложен в листинге 13.2. Большая часть этого листингаявляется вводом параметров задачи (шагов, начальных и граничных условий), и только в последней его строке определяется функция пользователя, вычисляющая сеточную функцию на каждом временном слое (при помощи встроенной функции решения системы линейных уравненийlsoive). В нескольких предыдущих строках листинга (после расчета коэффициента Куранта сои) формируется матрица системы уравнений А, которая записывается в подходящем для Mathcad виде, как это сделано в листинге 13.2. Как несложно убедиться, столбец правых частей разностныхуравнений выражается вычисленными значениями сеточной функции спредыдущего слоя.Результаты расчетов по неявной схеме показаны на рис. 13.12 и, как видно,они дают примерно те же результаты, что и в случае применения явнойсхемы (см.
рис. 13.7). Обратите внимание, что решение устойчиво при любых значениях коэффициента Куранта (в том числе, и больших 1), поскольку, как следует из соответствующих положений теории численных методов,неявная схема является безусловно устойчивой.Глава 13. Уравнения в частных производных333ф(х, u) :=0D:=lBorder (t) :=0Init (х) :=Ф(х-0.45) - Ф ( х - 0 . 5 5 )т:= 0.005М;=201Д:=—М= 0.05m:=0.. Mum:= Init(m2-D-TСои:=Сои = 4СоиСоит,т+1:~ААм,м:=10 , 0 :=V ( n ) :=uif n = Оl s o l v e ( A , - V { n - 1)otherwiseАлгоритм прогонкиПриведем в данном разделе описание чрезвычайно популярного алгоритмареализации неявных разностных схем, который называется методом прогонки.
Этот алгоритм имел историческое значение для становления технологийрасчетов уравнений в частных производных, и мы просто не можем не упомянуть о нем в этой книге.(ПримечаниеJСразу оговоримся, что его применение для решения уравнений в частных производных в среде Mathcad может быть оправдано, только если Вы работаете сочень частыми сетками, которые приводят к системам разностных уравненийбольшой размерности и, соответственно, очень долгому времени вычислений.Основным вычислительным ядром программы, реализующей на Mathcadнеявную разностную схему, было решение (на каждом временном слое)системы линейных алгебраических уравнений, задаваемых матрицей А. Заметим, что эта матрица, как говорят, имеет диагональное преобладание, аЧасть III.
Численные методы334точнее, является трехдиагональной (рис. 13.13). Все ее элементы, кроме элементов на главной диагонали и двух соседних диагоналях, равны нулю. Сточки зрения оптимизации быстродействия алгоритма, применение встроенной функции isoive является весьма расточительным, поскольку основной объем арифметических операций, выполняемых компьютером (а он составляет, как нетрудно убедиться величину порядка мг), сводится кнепроизводительному перемножению нулей.!о. г0-4i(3,<5Рис.
13.12. Решение линейного уравнения теплопроводностипри помощи неявной схемы на первом слое по времени (листинг 13.2)Для отыскания решения линейных систем алгебраических уравнений имеется чрезвычайно эффективный алгоритм, называемый прогонкой, которыйпозволяет снизить число арифметических операций на целый порядок, т. е.до значения порядка м. Это означает, что при использовании пространственных сеток с юоо узлами выигрыш во времени вычислений составит величину порядка ю 3 ! Реализация данного алгоритма приведена в листинге 13.3, который является продолжением листинга 13.2, используяопределенные в нем коэффициенты матрицы А, а также начальное условие.Программа листинга 13.3 осуществляет пересчет одного шага по времени,т.
е. заменяет содержимое столбца и с предыдущего временного слоя вычисленными значениями неизвестной функции со следующего слоя. Первыепять строк листинга 13.3 представляют так называемый обратный ход прогонки, а последние две строки — ее прямой ход. Заинтересовавшемуся читателю предлагается самому оформить представленный алгоритм прогонки ввиде программы решения разностных уравнений для вычисления произвольного временного слоя по примеру листингов 13.1 и 13.2. Заметим, чтоописание этого знаменитого алгоритма можно отыскать практически в любом современном учебнике по численным методам.Глава 13. Уравнения в частных производных33511to05•7Of 0UjOjDj•? 050fl00.51 4-2 1 0 5о5ьOf 0.50 00о_0_51Рис. 13.13. Матрица системы линейных разностных уравненийдля неявной схемы (листинг 13.2 для м=10)Листинг 13.3.
Алгоритм прогонки (продолжение листинга 13.2)ш.Рм-1 -СЕМ-1 := 0т:=М-1..У т :Ат,т +Ат , пн-1'iACtm-1 : = Y m m , m - lРт-1 : = 7 т - ( л ,т:=0..иМ-1пн-1:=ат-ит+13.2.3. О возможности решениямногомерных уравненийВсе, что было сказано до сих пор, касалось исключительно способов решения одномерных (в смысле пространственных координат) уравнений. И алгоритмы разностных схем, и встроенные функции, включая появившиеся в11-й версии (см. следующий разд.), относились к уравнениям, зависящим отодной пространственной координаты.Можно ли при помощи Mathcad решать двумерные или трехмерные (пространственные) уравнения? С точки зрения программирования пользователем численных алгоритмов типа метода сеток, принципиальных ограничений нет.
Разумеется, если сначала аккуратно выписать разностную схему336Часть HI. Численные методысоответствующего многомерного дифференциального уравнения, то вполневозможно запрограммировать ее при помощи описанных нами средств. Самым главным противодействием будет существенное увеличение временирасчетов. Простая оценка необходимого количества операций показывает,что ввод зависимости уравнения от второй пространственной координатымногократно увеличивает число разностных уравнений, которые должнырешаться при реализации каждого шага по времени. К примеру, если используется пространственная сетка из 100 узлов по каждой координате, товместо 102 разностных уравнений на каждом шаге придется решать уже 104уравнений, т. е. объем вычислений сразу же возрастает в 100 раз.
Вообщеговоря, пакет Mathcad не является экономичной средой вычислений, и бороться с их сильно возрастающим объемом пользователю следует еще наэтапе разработки алгоритма. Хорошим примером такой борьбы может служить применение специфических алгоритмов, типа метода прогонки (см.разд. «Алгоритм прогонки» этой главы).Приведем некоторые дополнительные замечания, связанные с возможностью осуществить редукцию громоздких (в смысле организации вычислений) двумерных задач к более простым. Рассмотрим ради примера двумерное уравнение теплопроводности (1) без источника с нулевыми граничнымиусловиями и некоторым начальным двумерным распределением температуры по расчетной поверхности:ди(х, у, t) _ Jd2u(x, у, t)Э;и(х, у, t)Произведем дискретизацию данного уравнения по временной координате,заменяя первую производную ее разностным аналогом и несколько перегруппировывая слагаемые и множители:ЭЧ.1Эх2Э+Ч.1ду2u=ui*i - i •DTКак Вы видите, мы используем неявную разностную схему, заранее заботясьо том, чтобы разностная задача была более устойчивой.
Здесь и £ (х,у) — известная с предыдущего шага по времени функция двух пространственныхпеременных, a u i + 1 (x,y) — функция, подлежащая определению при реализации каждого шага по времени.Можно посмотреть на полученную задачу с несколько другой стороны — аименно как на дифференциальное уравнение относительно неизвестнойфункции двух переменных u i + 1 (x,y). Подчеркнем, что такое уравнение получается для каждого шага по времени, т. е.
для реализации всей разностной схемы требуется решить большое число таких уравнений.С предложенной точки зрения, на каждом временном шаге необходимо решить некоторое двумерное эллиптическое линейное уравнение, причем егограничные условия определяются граничными условиями исходной задачи.Глава 13. Уравнения в частных производных337Это уравнение очень похоже на уравнение Пуассона с той лишь разницей,что в его правую часть, описывающую источник, входит неизвестная функция (к счастью, линейно).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
