Кирьянов. Самоучитель MathCad 11. Кирьянов д в
 Скачать 10.75 Mb.
|
|
Рис. 11.18. Решение жесткой системы ОДУ химической кинетики методом Рунге-Кутты (листинг 11.12) Глава Обыкновенные дифференциальные уравнения Функции для решения жестких ОДУ Решение жестких систем дифференциальных уравнений можно осуществить только с помощью встроенных функций, аналогичных по действию семейству рассмотренных выше функций для обычных ОДУ — алгоритм для жестких систем ОДУ fb(yo, — алгоритм для жестких систем ОДУ fr — алгоритм Розенброка для жестких систем ОДУ; • — вектор начальных значений в точке to; • — начальная и конечная точки расчетам — число шагов численного метода — векторная функция F(t,y) размера задающая систему ОДУ - матричная функция размера составленная из вектора производных функции F(t,y) t (правый столбец) и ее якобиана (N левых столбцов). Как Вы можете заметить, для двух последних функций серьезным отличием от функций, решающих нежесткие системы, является добавление к стандартному набору параметров дополнительной матричной функции, задающей якобиан системы ОДУ. Решение выдается в виде матрицы, по форме идентичной аналогичным функциям решения нежестких задач Коши. Покажем действие этих алгоритмов на том же примере жесткой системы ОДУ химической кинетики (листинг 11.14). Обратите внимание, как следует представлять в данном случае якобиан, сравнив задание матричной функции в предпоследней строке листинга 11.14 с заданием якобиана из листинга 11.13. Лиртинг 11.14. Решение жесткой системы ОДУ химической кинетики ( t , y ) : = 1 0 296 Часть III. Численные методы , у) := = stiffb о 0 (yO - 0 0 . 0 . 1 1 50 - 1 0 , 20 , F , J) 0 У1 Расчеты показывают, что для получения того же результата (см. рис. оказалось достаточно в тысячу раз меньшего количества шагов численного алгоритма, чем для стандартного метода Примерно во столько же раз требуется меньше компьютерного времени на проведение расчетов. Стоит ли говорить, что, если Вы имеете дело с жесткими (в той или иной степени) системами, применение описанных специальных алгоритмов просто необходимо. Важно заметить, что до сих пор мы имели дело с примером не очень жесткой системы. Попробуйте вместо скоростей упомянутых химических реакций (см. разд. 11.5.1), и взять другие числа, например о.О5, ю 4 и ю, соответственно. Заметим, что такое соотношение скоростей часто встречается в прикладных задачах химической кинетики и определяет куда более жесткую систему ОДУ. Ее уже никак не удается решить стандартными методами, поскольку число шагов численного метода должно быть просто гигантским. А между тем, алгоритмы для жестких ОДУ справляются с этой задачей с легкостью (рис причем практически при тех же значениях шага, что были взяты в листинге 11.14. Обратите внимание, что порядки величины решений для концентраций различных веществ на рис. 11.19 различаются еще сильнее, чем в предыдущем (менее жестком) примере. Примечание Это еще раз доказывает, что одна и та же система ОДУ с различными коэффициентами может быть жесткой в разной степени. В частности, приведенный выше пример генератора Ван дер Поля с параметром — это уже пример жесткой задачи. В заключение приведем соответствующие встроенные функции, которые применяются для решения жестких систем ОДУ не на всем интервале, а только водной заданной точке F, k, to, — метод Булирша-Штера Розенброка Имена этих функций пишутся с маленькой буквы, а их действие и набор параметров полностью аналогичны рассмотренным нами ранее для функций, относящихся к решению в заданной точке нежестких систем Глава Обыкновенные дифференциальные уравнения 297 (см. разд. 11.3.2). Отличие заключается в специфике применяемого ма и необходимости задания матричной функции якобиана Решение более жесткой системы ОДУ химической кинетики методом ГЛАВА Краевые задачи В этой главе рассматриваются краевые задачи для обыкновенных дифференциальных уравнений (ОДУ. Средства Mathcad позволяют решать краевые задачи для систем ОДУ, в которых часть граничных условий поставлена в начальной точке интервала, а остальная часть — в его конечной точке (см. разд. 12.1). Также возможно решать уравнения с граничными условиями, поставленными в некоторой внутренней точке. Для решения краевых задач в Mathcad предусмотрены соответствующие встроенные функции, реализующие алгоритм пристрелки (см. разд. Краевые задачи во множестве практических приложений часто зависят от некоторого числового параметра. При этом решение существует не для всех его значений, а лишь для счетного их числа. Такие задачи называют задачами на собственные значения (см. разд. Несмотря на то, что, в отличие от задач Коши для ОДУ, вне предусмотрены встроенные функции для решения жестких краевых задач, сними все-таки можно справиться, применив программирование разностных схем (см. разд. 12.3). 12.1. Краевые задачи для ОДУ Постановка краевых задач для ОДУ отличается от задач Коши, рассмотренных в главе 11, тем, что граничные условия для них ставятся неводной начальной точке, а на обеих границах расчетного интервала. Если имеется система N обыкновенных дифференциальных уравнений первого порядка, то часть из N условий может быть поставлена на одной границе интервала, а оставшиеся условия — на противоположной границе. Примечание Дифференциальные уравнения высших порядков можно свести к эквивалентной системе ОДУ первого порядка (см. гл. 11). Часть III. Численные методы постановке краевых задач Чтобы лучше понять, что из себя представляют краевые задачи, рассмотрим их постановочную часть на конкретном физическом примере модели взаимодействия встречных световых пучков. Предположим, что надо определить распределение интенсивности оптического излучения в пространстве между источником (лазером) и зеркалом, заполненном некоторой средой (рис. 12.1). Будем считать, что от зеркала отражается часть падающего излучения (те. его коэффициент отражения равен R), а среда как поглощает излучение с коэффициентом ослабления ах, таки рассеивает его. Причем коэффициент рассеяния назад равен В этом случае закон изменения интенсивности излучения, распространяющегося вправо, и интенсивности излучения влево определяется системой двух ОДУ первого порядка • + • . (1) dx = ах) • - r(x) Для правильной постановки задачи требуется, помимо уравнений, задать такое же количество граничных условий. Одно из них будет выражать известную интенсивность излучения ю, падающего с левой границы хо, а второе — закон отражения на его правой границе R (2) лазер 0 зеркало — свет У 1 X 12.1. Модель для постановки краевой задачи Полученную задачу называют краевой (boundary value problem), поскольку условия поставлены не на одной, а на обеих границах интервала Ив связи с этим, их не решить методами предыдущей главы, предназначенными для задач с начальными условиями. Далее для показа возможностей будем использовать этот пример с И конкретным видом и описывающим случай изотропного (независящего от координаты х) рассеяния Глава 12. Краевые Примечание Модель рис. 12.1 привела к краевой задаче для системы линейных ОДУ. Она имеет аналитическое решение в виде комбинации экспонент. Более сложные, нелинейные задачи, возможно решить только численно. Нетрудно сообразить, что модель станет нелинейной, если сделать коэффициенты ослабления и рассеяния зависящими от интенсивности излучения. Физически это будет соответствовать изменению оптических свойств среды под действием мощного излуче- ния. Примечание Модель встречных световых пучков привела нас к системе уравнений (1), в которые входят производные только по одной переменной х. Если бы мы стали рассматривать более сложные эффекты рассеяния в стороны (а не только впереди назад, тов уравнениях появились бы частные производные по другим пространственным переменным у и В этом случае получилась бы краевая задача для уравнений в частных производных, решение которой во много раз сложнее ОДУ. Алгоритм стрельбы Для решения краевых задач в Mathcad реализован наиболее популярный алгоритм, называемый методом стрельбы или пристрелки (shooting Он, по сути, сводит решение краевой задачи к решению серии задач Коши с различными начальными условиями. Рассмотрим здесь его основной принцип на примере модели (риса встроенные функции, реализующие этот алгоритм, приведем в следующем разделе. Суть метода стрельбы заключается в пробном задании недостающих граничных условий на левой границе интервала и решении затем полученной задачи Коши хорошо известными методами (см. гл. 11). В нашем примере не хватает начального условия для поэтому сначала зададим ему произвольное значение, например Конечно, такой выбор не совсем случаен, поскольку из физических соображений ясно, что, во- первых, интенсивность излучения — величина заведомо положительная, и, во-вторых, отраженное излучение должно быть намного меньше падающего. Решение задачи Коши с помощью функции приведено в листин- Листинг Решение пробной задачи Коши для модели \ ( ха ха х } + D ( х , ух Часть Численные методы := у 10 10 S r k f i x e d . 0 , 1 , График полученных решений показан на рис. 12.2 (слева. Из него видно, что взятое наугад второе начальное условие не обеспечило выполнение ничного условия при И понятно, что для лучшего выполнения этого граничного условия следует взять большее значение Возьмем, например и вновь решим задачу Результат показан на том же 12.2 (в центре. Граничное условие выполняется с лучшей точностью, но опять-таки оказалось недостаточным. Для еще одного значения получается решение, показанное на рис. 12.2 (справа. Из сравнения двух правых графиков легко заключить, что недостающее начальное условие больше 15, но меньше 20. Продолжая подобным образом "пристрелку" по недостающему начальному условию, возможно отыскать правильное решение краевой задачи. В этом и состоит принцип алгоритма стрельбы. Выбирая пробные начальные условия (проводя пристрелку) и решая соответствующую серию задач Коши, можно найти то решение системы ОДУ, которое (с заданной точностью) удовлетворит фаничному условию (или, в общем случае, условиям) на другой фанице расчетного интервала. (г) 10 0 . 5 0 . 5 . (о 0 0 . 5 12.2. Иллюстрация метода стрельбы (листинг Конечно, описанный алгоритм несложно запрограммировать самому, мив его как решение системы заданных алгоритмически уравнений, выражающих граничные условия на второй границе, относительно неизвестных пристрелочных начальных условий. Но делать этого нет необходимости, поскольку он оформлен в Mathcad в виде встроенных функций Глава 12. Краевые задачи 303 12.1.3. Решение двухточечных краевых задач Решение краевых задач для систем ОДУ методом стрельбы в Mathcad достигается применением двух встроенных функций. Первая предназначена для двухточечных задач с краевыми условиями, заданными на концах интервала load, score) — поиск вектора недостающих начальных условий для двухточечной краевой задачи для системы N ОДУ — вектор размера присваивающий недостающим начальным условиям (на левой границе интервала) начальные значения; • хО — левая граница расчетного интервала — правая граница расчетного интервала — векторная функция размера левых граничных условий, причем недостающие начальные условия именуются соответствующими компонентами векторного аргумента z; • у) — векторная функция размера выражающая L правых граничных условий для векторной функции у в точке — векторная функция, описывающая систему ОДУ, размера и двух аргументов — скалярного хи векторного у. При этому это неизвестная векторная функция аргументах того же размера Внимание! Решение краевых задач в Mathcad, в отличие от большинства остальных операций, реализовано не совсем очевидным образом. В частности, помните, что число элементов векторов D и l o a d равно количеству уравнений N, а векторов z, score и результата действия функции — количеству правых граничных условий L. Соответственно, левых граничных условий в задаче должно Как видно, функция sbval предназначена не для поиска собственно решения, те. неизвестных функций а для определения недостающих начальных условий впервой точке интервала, те Чтобы вычислить на всем интервале, требуется дополнительно решить задачу Коши. Разберем особенности использования функции sbval на конкретном примере (листинг 12.2), описанном выше (см. разд Краевая задача состоит из системы двух уравнений (N=2), ОДНОГО левого И ОДНОГО правого граничного условия Листинг 12.2. Решение краевой задачи ' ' I . 1 ах Часть III. Численные методы ( х ) ( х ) D ( х , ух х го := Ю 100 l o a d , z) s c o r e ( , у) :- У := s b v a l ( z , 0 , 1 , D , l o a d , s c o r e ( 18.555 > S , 1 , 10 , Первые три строки листинга задают необходимые параметры задачи и саму систему ОДУ. В четвертой строке определяется вектор z. Поскольку правое граничное условие всего одно, то недостающее начальное условие тоже одно, соответственно, и вектор z имеет только один элемент Ему необходимо присвоить начальное значение (мы приняли как в листинге, чтобы запустить алгоритм стрельбы (см. разд 12.1.2). Примечание Начальное значение фактически является параметром численного метода и поэтому может сильно повлиять на решение краевой задачи. В следующей строке листинга векторной функции присваиваются левые граничные условия. Эта функция аналогична векторной переменной, определяющей начальные условия для встроенных функций, решающих задачи Отличие заключается в записи недостающих условий. Вместо конкретных чисел на соответствующих местах пишутся имена искомых элементов вектора z. В нашем случае вместо второго начального условия стоит аргумент функции load. Первый аргумент функции load — это точка, в которой ставится левое граничное условие. Ее конкретное значение определяется непосредственно в списке аргументов функции Следующая строка листинга определяет правое граничное условие, для введения которого используется функция Оно записывается точно также, как система уравнений в функции D. Аргумент функции score аналогичен функции load и нужен для тех случаев, когда граничное условие явно зависит от координаты х. Вектор score должен состоять из такого же числа элементов, что и вектор Реализованный в функции sbval алгоритм стрельбы ищет недостающие начальные условия таким образом, чтобы решение полученной задачи Коши делало функцию как можно ближе к нулю. Как видно из листинга, результат применения sbval для интервала присваивается век Глава 12. Краевые задачи 305 торной переменной и. Этот вектор похож на вектор z, только в нем содержатся искомые начальные условия вместо приближенных начальных значений, заданных в z. Вектор содержит, как и z, всего один элемент С его помощью можно определить решение краевой задачи ух) (последняя строка листинга. Тем самым, функция сводит решение краевых задач к задачам График решения краевой задачи показан на рис. 12.3, 12.3. Решение краевой задачи для (листинг Рис. 12.4. Решение краевой задачи для На рис показано решение той же самой краевой задачи, нос другим правым граничным условием, соответствующим Те. без зеркала на правой границе. В этом случае слабый обратный пучок света образуется исключительно за счет обратного рассеяния излучения от лазера. Конечно, многие из читателей уже обратили внимание, что реальная физическая среда не может создавать такого большого рассеяния назад. Иными словами, более реальны значения Однако, когда коэффициенты в системе ОДУ при разных очень сильно (на порядки) различаются, система ОДУ становится жесткой, и функция sbval не может найти решения, выдавая вместо него сообщение об not find a solution"). Внимание! Метод стрельбы не годится для решения жестких краевых задач. Поэтому алгоритмы решения жестких ОДУ в Mathcad приходится программировать самому (см. разд. 12.3). 12.1.4. Решение краевых задач с дополнительным условием в промежуточной точке Иногда дифференциальные уравнения определяются с граничными условиями не только на концах интервала, но и с дополнительным условием в 306 Часть III. Численные методы некоторой промежуточной точке расчетного интервала. Чаще всего такие задачи содержат данные о негладких в некоторой внутренней точке интервала решениях. Для них имеется встроенная функция bvaifit, также реализующая алгоритм стрельбы — ПОИСК вектора достающих граничных условий для краевой задачи с дополнительным условием в промежуточной точке для системы ОДУ — вектор, присваивающий недостающим начальным условиям на левой границе интервала начальные значения вектор того же размера, присваивающий недостающим начальным условиям на правой границе интервала начальные значения — левая граница расчетного интервала — правая граница расчетного интервала — точка внутри интервала векторная функция, описывающая систему N ОДУ, размера и двух аргументов — скалярного хи векторного у. При этому это неизвестная векторная функция аргументах того же размера — векторная функция размера левых граничных условий, причем недостающие начальные условия поименовываются соответствующими компонентами векторного аргумента z; • — векторная функция размера правых граничных условий, причем недостающие начальные условия поименовываются соответствующими компонентами векторного аргумента — векторная функция размера выражающая внутреннее условие для векторной функции у в точке Обычно функция bvaifit применяется для задач, в которых производная решения имеет разрыв во внутренней точке xf. Некоторые из таких задач не удается решить обычным методом пристрелки, поэтому приходится вести пристрелку сразу из обеих граничных точек. В этом случае дополнительное внутреннее условие в точке xf является просто условием сшивки в ней левого и правого решений. Поэтому для данной постановки задачи оно записывается в форме у) :=у. Рассмотрим действие функции bvaifit на знакомом примере модели взаимодействия пучков света (см. рис предположив, что в промежутке между и находится другая, оптически более плотная среда с другим коэффициентом ослабления излучения ах) з. Соответствующая краевая задача решена в листинге 12.3, причем разрывный показатель ослабления определяется в его второй строке Глава 12. Краевые задачи 307 |