Кирьянов Д. - MathCad 11 (1077323), страница 48
Текст из файла (страница 48)
При меньшем числе шагов численному алгоритму неудается найти решение. В процессе работы алгоритма оно расходится, иMathcad вместо результата выдает ошибку о превышении предельно большого числа.Еще один факт, на который стоит обратить внимание, — это различие впорядке величины получающегося решения.
Как видно из рис. 11.18, концентрация первого реагента у : существенно (в тысячи раз) превышает концентрацию остальных. Это свойство также очень характерно для жесткихсистем.ПримечаниеJВ принципе, можно было бы снизить жесткость системы "вручную", применяямасштабирование. Для этого нужно искусственно уменьшить искомую функциюy l , к примеру, в тысячу раз, разделив все слагаемые в системе ОДУ, содержащие y i , на 1000.
После масштабирования для решения полученной системыметодом Рунге-Кутты будет достаточно взять всего М-20 шагов.Рис. 11.18. Решение жесткой системы ОДУ химической кинетики методомРунге-Кутты (листинг 11.12)Глава 11. Обыкновенные дифференциальные уравнения29511.5.2. Функции для решения жестких ОДУРешение жестких систем дифференциальных уравнений можно осуществитьтолько с помощью встроенных функций, аналогичных по действию семейству рассмотренных выше функций для обычных ОДУ.П Radau(yo, to,ti,M,F) — алгоритм RADAUS для жестких систем ОДУ;Оs t i f fb(yo, to, ti,M,F, j ) — алгоритм Булирша-Штера для жестких систем ОДУ;Пs t i f fr (yo, to, ti,M,F, j ) — алгоритм Розенброка для жестких системОДУ;• уо — вектор начальных значений в точке to;•t o , t i — начальная и конечная точки расчета;• м — число шагов численного метода;• F — векторная функция F(t,y) размера IXN, задающая систему ОДУ;•j - матричная функция j ( t , y ) размера (N+DXN, составленная извектора производных функции F(t,y) no t (правый столбец) и ееякобиана (N левых столбцов).Как Вы можете заметить, для двух последних функций серьезным отличиемот функций, решающих нежесткие системы, является добавление к стандартному набору параметров дополнительной матричной функции, задающей якобиан системы ОДУ.
Решение выдается в виде матрицы, по формеидентичной аналогичным функциям решения нежестких задач Коши.Покажем действие этих алгоритмов на том же примере жесткой системыОДУ химической кинетики (листинг 11.14). Обратите внимание, как следуетпредставлять в данном случае якобиан, сравнив задание матричной функции в предпоследней строке листинга 11.14 с заданием якобиана из листинга 11.13.Лиртинг 11.14.
Решение жесткой системы ОДУ химической кинетикиуО:=- 0 . 1 - у 0 + 10F ( t , y ) :='У1-У20.1-Уо- Ю2-у1-уа " 10J103-У1296Часть III. Численные методыf оJt , у) :=-0100. 100D = stiffb (yO. о.ioio2-y2-102-y2-10io332-ioУ12•Yl050 , 20 , F , J )Расчеты показывают, что для получения того же результата (см. рис. 11.18)оказалось достаточно в тысячу раз меньшего количества шагов численногоалгоритма, чем для стандартного метода Рунге-Кутты! Примерно во столькоже раз требуется меньше компьютерного времени на проведение расчетов.Стоит ли говорить, что, если Вы имеете дело с жесткими (в той или инойстепени) системами, применение описанных специальных алгоритмов просто необходимо.Важно заметить, что до сих пор мы имели дело с примером не очень жесткой системы.
Попробуйте вместо скоростей упомянутых химических реакций (см. разд. 11.5.1), o.i, ю 3 и ю 2 взять другие числа, например о.О5, ю4и ю 7 , соответственно. Заметим, что такое соотношение скоростей частовстречается в прикладных задачах химической кинетики и определяет кудаболее жесткую систему ОДУ. Ее уже никак не удается решить стандартнымиметодами, поскольку число шагов численного метода должно быть простогигантским. А между тем, алгоритмы для жестких ОДУ справляются с этойзадачей с легкостью (рис. 11.19), причем практически при тех же значенияхшага, что были взяты в листинге 11.14.
Обратите внимание, что порядкивеличины решений для концентраций различных веществ на рис. 11.19 различаются еще сильнее, чем в предыдущем (менее жестком) примере.ПримечаниеЭто еще раз доказывает, что одна и та же система ОДУ с различными коэффициентами может быть жесткой в разной степени. В частности, приведенныйвыше пример генератора Ван дер Поля с параметром |л=5000 — это уже пример жесткой задачи.В заключение приведем соответствующие встроенные функции, которыеприменяются для решения жестких систем ОДУ не на всем интервале, атолько в одной заданной точке t i .• radau(y0, t o , t l , асе, F, k, s) — метод RADAUSПstiffb(yo, to, t i , a c c , F , j , k , s ) — метод Булирша-ШтераПs t i f fr (yo, to, t i , a c c , F , j , k , s ) — метод РозенброкаИмена этих функций пишутся с маленькой буквы, а их действие и наборпараметров полностью аналогичны рассмотренным нами ранее для функций, относящихся к решению в заданной точке нежестких системГлава 11.
Обыкновенные дифференциальные уравнения297(см. разд. 11.3.2). Отличие заключается в специфике применяемого алгоритма и необходимости задания матричной функции якобиана j ( c , y ) .Рис. 11.19. Решение более жесткой системы ОДУ химической кинетикиметодом РозенброкаГЛАВА 12Краевые задачиВ этой главе рассматриваются краевые задачи для обыкновенных дифференциальных уравнений (ОДУ).
Средства Mathcad позволяют решать краевые задачи для систем ОДУ, в которых часть граничных условий поставлена в начальной точке интервала, а остальная часть — в его конечной точке(см. разд. 12.1). Также возможно решать уравнения с граничными условиями, поставленными в некоторой внутренней точке. Для решения краевыхзадач в Mathcad предусмотрены соответствующие встроенные функции,реализующие алгоритм пристрелки (см.
разд. 12.1.2).Краевые задачи во множестве практических приложений часто зависят отнекоторого числового параметра. При этом решение существует не для всехего значений, а лишь для счетного их числа. Такие задачи называют задачами на собственные значения (см. разд. 12.2).Несмотря на то, что, в отличие от задач Коши для ОДУ, в Mathcad не предусмотрены встроенные функции для решения жестких краевых задач, с нимивсе-таки можно справиться, применив программирование разностных схем(см. разд. 12.3).12.1.
Краевые задачи для ОДУПостановка краевых задач для ОДУ отличается от задач Коши, рассмотренных в главе 11, тем, что граничные условия для них ставятся не в одной начальной точке, а на обеих границах расчетного интервала. Если имеетсясистема N обыкновенных дифференциальных уравнений первого порядка,то часть из N условий может быть поставлена на одной границе интервала, аоставшиеся условия — на противоположной границе.ПримечаниеjДифференциальные уравнения высших порядков можно свести к эквивалентной системе ОДУ первого порядка (см. гл.
11).Часть III. Численные методы30012.1.1.0 постановке краевых задачЧтобы лучше понять, что из себя представляют краевые задачи, рассмотримих постановочную часть на конкретном физическом примере модели взаимодействия встречных световых пучков. Предположим, что надо определитьраспределение интенсивности оптического излучения в пространстве междуисточником (лазером) и зеркалом, заполненном некоторой средой(рис. 12.1).
Будем считать, что от зеркала отражается R-Я часть падающегоизлучения (т. е. его коэффициент отражения равен R), а среда как поглощает излучение с коэффициентом ослабления а(х), так и рассеивает его. Причем коэффициент рассеяния назад равен г(х). В этом случае закон изменения интенсивности у о (х) излучения, распространяющегося вправо, иинтенсивности yi(x) излучения влево определяется системой двух ОДУпервого порядка:0= ~а(х) • у о (х> + г(х> • У 1 {х)dYl(x)= а(х) •dx.(1)- r(x) • y o (x)Для правильной постановки задачи требуется, помимо уравнений, задатьтакое же количество граничных условий. Одно из них будет выражать известную интенсивность излучения ю, падающего с левой границы х=о, авторое — закон отражения на его правой границе x=i:у о (О) = 10(2)= R-зеркало^лазер————г-свет//1УУ1101XРис.
1 2 . 1 . Модель для постановки краевой задачиПолученную задачу называют краевой (boundary value problem), поскольку условия поставлены не на одной, а на обеих границах интервала {0,1).И, в связи с этим, их не решить методами предыдущей главы, предназначенными для задач с начальными условиями. Далее для показа возможностей Mathcad будем использовать этот пример с R=I И конкретным видомa(x)=const=i и г(x)=const=o.i, описывающим случай изотропного (не зависящего от координаты х) рассеяния.Глава12.Краевые задачи301ПримечаниеМодель рис.
12.1 привела к краевой задаче для системы линейных ОДУ. Онаимеет аналитическое решение в виде комбинации экспонент. Более сложные,нелинейные задачи, возможно решить только численно. Нетрудно сообразить,что модель станет нелинейной, если сделать коэффициенты ослабления и рассеяния зависящими от интенсивности излучения. Физически это будет соответствовать изменению оптических свойств среды под действием мощного излучения.ПримечаниеМодель встречных световых пучков привела нас к системе уравнений (1), в которые входят производные только по одной переменной х.
Если бы мы сталирассматривать более сложные эффекты рассеяния в стороны (а не только вперед и назад), то в уравнениях появились бы частные производные по другимпространственным переменным у и ъ. В этом случае получилась бы краеваязадача для уравнений в частных производных, решение которой во много разсложнее ОДУ.12.1.2. Алгоритм стрельбыДля решения краевых задач в Mathcad реализован наиболее популярныйалгоритм, называемый методом стрельбы или пристрелки (shooting method).Он, по сути, сводит решение краевой задачи к решению серии задач Кошис различными начальными условиями. Рассмотрим здесь его основнойпринцип на примере модели (рис. 12.1), а встроенные функции, реализующие этот алгоритм, приведем в следующем разделе.Суть метода стрельбы заключается в пробном задании недостающих граничных условий на левой границе интервала и решении затем полученнойзадачи Коши хорошо известными методами (см.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
