Кирьянов. Самоучитель MathCad 11. Кирьянов д в
 Скачать 10.75 Mb.
|
|
Листинг 12.3. Краевая задача условием в промежуточной точке := 100 := (х) := .1 1 D(x, уха (х + хау у ( zl , z2 , 0 , 5 , loadl , , score) ( Система уравнений и левое краевое условие вводится также, как ив предыдущем листинге для функции Обратите внимание, что таким же образом записано и правое краевое условие. Для того чтобы ввести условие отражения на правой границе, пришлось определить еще один неизвестный пристрелочный параметр Строка листинга, в которой определена функция score, задает условие стрельбы — сшивку двух решений в точке В самой последней строке листинга выдан ответ — определенные численным методом значения обоих пристрелочных параметров, которые объединены в вектор (мы применили в предпоследней строке операцию транспонирования, чтобы результат получился в форме вектора, а не мат- рицы-строки). Для корректного построения фафика решения лучше составить его из двух частей — решения задачи на интервале и другой задачи ши для интервала Реализация этого способа приведена в листинге, который является продолжением листинга 12.3. В последней строке листинга 12.4 выведено значение второй искомой функции на правой границе интервала. Всегда полезно проконтролировать, что оно совпадает с соответствующим пристрелочным параметром в последней строке листинга 12.3). Часть HI. Численные методы Листинг 12.4. Решение краевой задачи (продолжение листинга := 10 V SO rkfixed := , 0 , 0.5 , , D , 0.5 , 1 , M, Решение краевой задачи приведено на рис. 12.5. С физической точки зрения естественно, что интенсивность света уменьшается быстрее по мере распространения в более плотной среде в правой половине расчетного интервала. В средней точке как и ожидалось, производные обоих решений имеют разрыв. Примечание Еще один пример решения краевых задач с разрывными коэффициентами ОДУ приведен в справочной системе Рис. 12.5. Решение краевой задачи с разрывом в x f = 0 . листинги Ради справедливости необходимо заметить, что разобранную краевую задачу легко решить и с помощью функции заменив в листинге 12.2 зависимость ах) на третью строку листинга 12.3. В этом случае листинг даст в точности тот же ответ, что показам на рис. 12.5. Однако в определенных (в том числе из соображений быстроты расчетов) удобнее использовать функцию те. вести пристрелку с обеих границ интер- вата Глава Краевые задачи 309 Совет Если Вы имеете дело с подобными уравнениями, попробуйте сначала решить их как обычную краевую задачу с помощью более надежной и легкой в применении функции Задачи на собственные значения для ОДУ Задачи на собственные значения — это краевые задачи для системы ОДУ, в которой правые части зависят от одного или нескольких параметров Значения этих параметров неизвестны, а решение краевой задачи существует только при определенных которые называются собственными значениями (eigenvalues) задачи. Решения, соответствующие этим называют собственными функциями (eigenfunctions) задачи. Правильная постановка таких задач требует формулировки количества граничных условий, равного сумме числа уравнений и числа собственных значений. Физическими примерами задач на собственные значения являются, например, уравнение колебаний струны, уравнение Шредингера в квантовой механике, уравнения волн в резонаторах и многие другие. С вычислительной точки зрения, задачи на собственные значения очень похожи на рассмотренные выше краевые задачи. В частности, для многих из них также применим метод стрельбы (см. разд. 12.1.2). Отличие заключается в пристрелке не только по недостающим левым граничным условиям, но еще и по искомым собственным значениям. В Mathcad для решения задач на собственные значения используются те же функции sbval и В их первый аргумент, те. вектор, присваивающий начальные значения недостающим начальным условиям, следует включить и начальное приближение для собственного значения. Рассмотрим методику решения на конкретном примере определения собственных упругих колебаний струны. Профиль колебаний струны ух) описывается линейным дифференциальным уравнением второго порядка dx dx где и q(x) жесткость и плотность, которые, вообще говоря, могут меняться вдоль струны. Если струна закреплена на обоих концах, то граничные условия задаются в виде Сформулированная задача является частным случаем задачи Поскольку решается система двух ОДУ, содержащая одно собственное значение X, то по идее задача требует задания трех (2+1) условий. Однако, как легко убедиться, уравнение колебаний струны — линейное и однородное, поэтому в любом случае решение ух) будет определено с точностью до множителя. Это Часть III. Численные методы что производную решения можно задать произвольно, например у что и будет третьим условием. Тогда краевую задачу можно решать как задачу а пристрелку вести только по одному параметру — венному значению. Процедура поиска первого собственного значения представлена в листинге Листинг 12.5. Решение задачи о собственных колебаниях р ( х ) := 1 ( х ) := 1 р ' ( х ) О := О Ь := 1 := 0 . 5 D ( х , у ) := l o a d ( a , Я ) : - У1 Р ( х ; ( • ( р ' ( х ) + •; s c o r e (Ь , у уо Л :- , a b , D , l o a d , s c o r e Л ( 9.87 = В первых двух строках листинга определяются функции, входящие в задачу, в том числе и границы расчетного интервала В третьей строке дается начальное приближение к собственному значению в четвертой вводится система ОДУ. Обратите внимание, что она состоит не из двух, а из трех уравнений. Первые два из них определяют эквивалентную (систему ОДУ первого порядка, а третье необходимо для задания собственного значения в виде еще одного компонента искомого вектора у. Поскольку, по определению, собственное значение постоянно при всех х, то его производная должна быть приравнена нулю, что отражено в последнем уравнении. Важно также, что во втором из уравнений собственное значение записано как поскольку является одним из неизвестных. В следующих двух строках листинга задается левое граничное условие, включающее и недостающее условие на собственное значение для третьего уравнения, и правое граничное условие В предпоследней строке листинга обычным образом применяется функция sbval, а в последней выводится результат ее работы вместе с известным аналитически собственным значением Как легко убедиться, мы нашли первое собственное значе- Глава 12. Краевые задачи 311 ние для а чтобы найти другие собственные значения, необходимо задать другие начальные приближения к ним (в третьей строке листинга Например, выбор приводит ко второму собственному значению а — к третьему Чтобы построить график соответствующей собственной функции, надо добавить в листинг строку, программирующую решение задачи например, такую Полученные кривые показаны на рис. 12.6 в виде коллажа трех графиков, рассчитанных для трех собственных значений. Примечание Примеры решения нескольких задач на собственные значения можно найти в разделе Mathcad Рис. 12.6. Первые три собственные функции задачи колебаний струны (коллаж трех графиков Разностные схемы для ОДУ Многие краевые задачи не поддаются решению методом стрельбы. Однако в Mathcad других встроенных алгоритмов нет. Тем не менее, это не означает, что по-другому решать краевые задачи невозможно, ведь другие численные алгоритмы несложно запрограммировать самому пользователю. Рассмотрим возможную реализацию наглядного метода, называемого разностным, которым можно решать краевые задачи как для ОДУ, таки для дифференциальных уравнений в частных производных О разностном методе решения ОДУ Разберем идею разностного метода решения краевых задач на примере взаимодействия световых пучков (см. рис. 12.1), переобозначив в системе) интенсивность излучения вправо на у, а интенсивность излучения Зак 984 Часть III. Численные методы влево на у (просто в целях удобства, чтобы не писать индекс. Суть метода заключается в покрытии расчетного интервала сеткой из точек. Тем самым определяются шагов (рис Затем надо заменить дифференциальные уравнения исходной краевой задачи аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого шага. В нашем случае достаточно просто заменить первые производные из (12.1,1) их разностными аналогами (такой метод называется еще методом Эйлера): Д ( 1 ) • \/• — - V. 0 0 Y 1 У i +Д - 1 1 X 12.7. Сетка, покрывающая расчетный интервал Примечание Существует множество способов аппроксимации дифференциальных уравнений разностными. От выбора конкретного варианта зависит не только простота, быстрота и удобство вычислений, но и сама возможность получения правильного ответа. Получилась система (по числу шагов разностных линейных алгебраических уравнений с неизвестными и Для того чтобы она имела единственное решение, надо дополнить число уравнений до ЭТО МОЖНО сделать, записав в разностном виде оба граничных условия = I 0 , (Сформированная полная система алгебраических уравнений называется разностной схемой аппроксимирующей исходную краевую задачу. Обратите внимание, что правые части разностных уравнений системы (1) на каждом шаге записаны для левой границы шага. Такие разностные схемы называют явными т. к. все значения и находятся в левой части уравнений. Полученную явную разностную схему легко записать в матричной форме (где — неизвестный вектор, получающийся объединением векторов Y и у. Решив систему (3), мы получим решение краевой задачи Глава 12. Краевые задачи 313 Примечание На самом деле, все несколько сложнее, поскольку, вообще говоря, необходимо еще доказать, что, во-первых, разностная схема действительно аппроксимирует дифференциальные уравнения и, во-вторых, при разностное решение действительно сходится к дифференциальному. Процесс решения системы разностных уравнений называют также реализацией разностной схемы. Программа, которая решает рассматриваемую краевую задачу разностным методом, приведена в листинге 12.6. I Листинг Реализация явной разностной схемы ах := 1 , := ( i Д - i : = N + 1 ., := - - 1 , := 1 := Д • r ( i Д - 1 N-2 i := 1 .. 2 - N - l := 0 := A := , 0 , 2 - 1 , 1 , 2 := (A, B) = 13.453 = 3 1 . 7 Дадим минимальные комментарии, надеясь, что заинтересовавшийся читатель с карандашом в руках разберется в порядке индексов и соответствии матричных элементов, а возможно, составит и более удачную программу. В первой строке листинга определяются функции и константы, входящие в модель, во второй задается число точек сетки и ее равномерный шаг. Следующие две строки определяют матричные коэффициенты, аппроксимирующие уравнения для а пятая и шестая — для Седьмая и восьмая строки листинга задают, соответственно, левое и правое граничное условие, а строки с девятой по одиннадцатую правые части системы (3). В следующей строке завершается построение матрицы А вырезанием из нее левого нулевого столбца. В предпоследней строке листинга применена встроенная функция lsolve для решения системы (3), а в последней выведены рассчитанные ею неизвестные граничные значения. Графики решения при Часть HI. Численные методы ведены на рис. 12.8, причем первые N элементов итогового вектора есть вычисленное излучение вперед, а последние N элементов излучение назад. л- 50 а. ч 0 1 0.5 Д — 12.8. Решение краевой задачи разностным методом (листинг Как мы увидели, реализация в Mathcad разностных схем вполне возможна и не слишком трудоемка — предложенная программа состоит всего из двух десятков математических выражений. Конечно, для их написания требуется и время, и часто кропотливые расчеты, но, собственно, в этом и состоит работа математика. Кстати говоря, при небольшом числе шагов, расчеты по разностным схемам не требуют существенного времени (программа, приведенная в листинге 12.6, работает быстрее, чем метод стрельбы, встроенный в функцию Существуют, кроме того, весьма очевидные для многих читателей пути ускорения расчетов, связанные с применением более подходящих методов решения систем линейных уравнений с разреженной матрицей. Жесткие краевые задачи Один из случаев, когда применение разностных схем может быть очень лезным, связан с решением жестких краевых задач (подробнее о жестких ОДУ читайте в гл. 11). В частности, рассматриваемая задача о встречных световых пучках становится жесткой при увеличении коэффициента ослабления ах) в несколько десятков раз. Например, при попытке решить ее с ах) :=юо с помощью листинга 12.2 вместо ответа выдается сообщение об ошибке "Can't converge to a solution. Encountered too many integration steps" сходится к решению. Слишком много шагов интегрирования. Это и неудивительно, поскольку жесткие системы характерны тем, что требуют исключительно малого значения шага в стандартных алгоритмах. Для жестких задач неприменимы и явные разностные схемы, о которых рассказывалось в предыдущем разделе (см. разд. 12.3.1). Результат расчетов Глава 12. Краевые задачи 315 по программе листинга 12.6, например с (рис. 12.9), дает характерную для неустойчивых разностных схем "разболтку" — колебания нарастающей амплитуды, не имеющие ничего общего с реальным решением Неверное решение жесткой краевой задачи по неустойчивой явной разностной схеме Выходом из положения будет использование неявных разностных схем. Применительно к нашей задаче достаточно заменить правые части уравнений значениями не на левой, а на правой границе каждого шага + 1 + l) + + ( 4 Граничные условия, конечно, можно оставить в том же виде (2). Поскольку мы имеем дело с линейными дифференциальными уравнениями, то и схему) легко будет записать в виде матричного равенства (3), перегруппировывая соответствующим образом выражение (4) и приводя подобные слагаемые. Разумеется, полученная матрица А будет иной, нежели матрица А ДЛЯ явной схемы (1). Поэтому и решение (реализация неявной схемы) может отличаться от изображенного на рис. 12.9 результата расчетов по явной схеме. Программа, составленная для решения системы (4), приведена в листинге Листинг Реализация неявной разностной схемы для жесткой краевой задачи. ах Д := N - 1 i := 1 .. N - 1 Часть III. Численные методы , i := - 1 := -Д -Д - д : = N + 1 .. 2 : = , 1 N • 2 N i := 1 .. := О := ( A , 0 , 2 1 , 1 , 2 : = ( A , B ) . - 6 := a ( i - + 1 := Д 0 . 0 1 7 x 1 0 будем специально останавливаться на обсуждении листинга 12.7, поскольку он почтив точности повторяет предыдущий листинг. Отличие заключается лишь в формировании матрицы А другим способом, согласно неявной схеме. Решение, показанное на рис. 12.10, демонстрирует, что произошло чудо "разболтка" исчезла, а распределение интенсивностей стало физически предсказуемым. Обратите внимание, что (из-за взятого нами слишком большого коэффициента ослабления излучения) отраженный пучок света имеет очень маленькую интенсивность, и ее пришлось построить на графике с увеличением в тысячу раз ;= 0 . N - 1 12.10. Решение краевой задачи разностным методом 12.7) ГЛАВА 13 Дифференциальные уравнения в частных производных Дифференциальные уравнения в частных производных представляют собой одну из наиболее сложных и одновременно интересных задач вычислительной математики. Эти уравнения характеризуются тем, что для их решения не существует единого универсального алгоритма и большинство задач требует своего собственного особого подхода. Уравнениями в частных производных описывается множество разнообразных физических явлений, и сих помощью можно с успехом моделировать самые сложные явления и процессы (диффузия, гидродинамика, квантовая механика, экология и т. д.). Дифференциальные уравнения в частных производных требуют нахождения функции не одной, как для ОДУ, а нескольких переменных, например f(x,y) или f(x,t). Постановка задач (см. разд. 13.1) включает в себя само уравнение (или систему уравнений, содержащее производные неизвестной функции по различным переменным (частные производные, а также определенное количество краевых условий на границах расчетной области. Несмотря на то, что обладает довольно ограниченными возможностями по отношению к уравнениям в частных производных, в нем имеется несколько встроенных функций, количество и возможности которых увеличились в новой й версии (см. разд. 13.3). Решать уравнения в частных производных можно и путем непосредственного программирования пользовательских алгоритмов (см. разд. 13.2). Автор совершенно сознательно сначала рассматривает численные методы для решения уравнений в частных производных, а уже затем описывает предназначенные для этого встроенные функции, чтобы читатель ясно осознавал, каким образом Mathcad производит расчеты использование встроенных функций для решения уравнений в частных производных не всегда бывает успешными Часть Численные методы ответственность за верный выбор их параметров часто ложится на пользователя, которому необходимо четко представлять основные принципы функционирования численных алгоритмов, примененных во встроенных функ- циях. Рассмотрение уравнений в частных производных построено на анализе различных вариаций уравнения диффузии тепла, которое допускает наглядную физическую интерпретацию. Постановка задач. Классификация уравнений в частных производных Постановка задач для уравнений в частных производных включает определение самого уравнения (или системы нескольких уравнений, а также необходимого количества краевых условий (число и характер задания которых определяется спецификой уравнения. По своему названию уравнения должны содержать частные производные неизвестной функции (ИЛИ нескольких функций, если уравнений несколько) по различным аргументам, например пространственной переменной хи времени t. Соответственно, для решения задачи требуется вычислить функцию нескольких переменных, например в некоторой области определения аргументов ох Ио < t < т. Граничные условия определяются как заданные временные зависимости функции и, или производных этой функции на границах расчетной области о и а начальные — как заданная (х, о). Сами уравнения в частных производных (несколько условно) можно разделить натри основных типа параболические — содержащие первую производную по одной переменной и вторую подругой, причем все эти производные входят в уравнение с одинаковым знаком гиперболические - содержащие первую производную по одной переменной и вторую подругой, входящие в уравнение с разными знаками эллиптические — содержащие только вторые производные, причем одного знака. Некоторые более сложные уравнения нельзя однозначно подогнать под приведенную классификацию, тогда говорят о гибридных типах уравнений. Пример уравнение диффузии тепла На протяжении всей главы мы будем использовать в качестве примера очень наглядное и имеющее различные, от очевидных до самых неожиданных, решения уравнение теплопроводности Глава 13. Уравнения в частных производных 319 Двумерное динамическое уравнение Рассмотрим следующее параболическое уравнение в частных производных, зависящее от трех переменных — двух пространственных хи у, а также от времени t: t) _ Эх+ ф(х, у, t, Примечание Выражение в скобках в правой части уравнения (сумму вторых пространственных производных функции и часто, ради краткости, обозначают при помощи оператора Лапласа: Это уравнение называется двумерным уравнением теплопроводности или, по- другому уравнением диффузии тепла Оно описывает динамику распределения температуры на плоской поверхности (например, на металлической пластине) в зависимости от времени 13.1). Физический смысл коэффициента D, который, вообще говоря, может быть функцией как координат, таки самой температуры, заключается в задания скорости перетекания тепла от более нагретых областей в менее нагретые. Функция описывает приток тепла извне, те. источники тепла, которые также могут зависеть как и от пространственных координат (что задает локализацию источников, таки от времени и температуры и 13.1. модель, описываемая двумерным уравнением теплопроводности Для того чтобы правильно поставить краевую задачу для двумерного уравнения теплопроводности, следует определить следующие дополнительные условия: • граничные условия, те. динамику функции u(x,y,t) и / или ее производных на границах расчетной области; • начальное условие, те. функцию Часть Численные методы Примечание Если рассматривается не одно уравнение в частных производных, а система уравнений, то соответствующие начальные и граничные условия должны быть поставлены для каждой из неизвестных функций. Стационарное двумерное уравнение Частный случай уравнения теплопроводности определяет стационарную, т. е. независящую от времени задачу. Стационарное уравнение описывает физическую картину распределения температуры по пластине, не изменяющуюся стечением времени. Такая картина может возникнуть при условии, что стационарный источник тепла действует довольно продолжительное время, и переходные процессы, вызванные его включением, прекратились. Пример численного решения такого уравнения показан на рис. 13.2 в виде поверхности и(х,у). |