Кирьянов Д. - MathCad 11 (1077323), страница 45
Текст из файла (страница 45)
В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором столбце— первой компоненты и т. д.Чтобы построить график решения, надо отложить соответствующие компоненты матрицы решения по координатным осям: значения аргумента и<0> —вдоль оси х, а и*1' и и <2> — вдоль оси у (рис. 11.4).
Как известно, решенияобыкновенных дифференциальных уравнений часто удобнее изображать неЧасть III. Численные методы276в таком виде, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций. При этом аргумент входит в них лишь параметрически. В рассматриваемом случае двух ОДУ такойграфик — фазовый портрет системы — является кривой на фазовой плоскости и поэтому особенно нагляден. Он изображен на рис. 11.5 (слева), иможно заметить, что для его построения потребовалось лишь поменять метки осей на и*1" и и <2> , соответственно.ПримечаниеФазовый портрет типа, изображенного на рис.
11.5, имеет одну стационарнуюточку (аттрактор), на которую "накручивается" решение. В теории динамических систем аттрактор такого типа называется фокусом.Рис. 11.4. График решения системы ОДУ (11.2—1) (листинг 11.4)п-О.1-0.05г-О.1-0.050.050.1Рис. 1 1 . 5 . Фазовый портрет решения системы ОДУ при м=100 (слева)и М=200 (справа) (листинг 11.4)В общем случае, если система состоит из N ОДУ, то фазовое пространствоявляется N-мерным. При N>3 наглядность теряется, и для визуализации фазового портрета приходится строить его различные проекции.Глава 11.
Обыкновенные дифференциальные уравнения277На том же рис. 11.5, справа, показан для сравнения результат расчета фазового портрета с большим числом шагов. Видно, что в этом случае обеспечивается лучшая точность и, в результате, решение получается более гладким.Конечно, при этом увеличивается и время расчетов.Внимание!При решении систем ОДУ многие проблемы могут быть устранены при помощипростой попытки увеличить число шагов численного метода. В частности, сделайте так при неожиданном возникновении ошибки "Found a number with amagnitude greater than 10Л307" (Найдено число, превышающее значение10307).
Данная ошибка может означать не то, что решение в действительностирасходится, а просто недостаток шагов для корректной работы численного алгоритма.В заключение следует сказать несколько слов об особенностях различныхчисленных методов. Все они основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации, получаются алгоритмы различной точности и быстродействия. В Mathcad использован наиболее популярный алгоритм РунгеКутты четвертого порядка, описанный в большинстве книг по методам вычислений.
Он обеспечивает малую погрешность для широкого класса системОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо rkfixed другие функции, благо сделать это очень просто,благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе.Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо, либо существуют участки медленных и быстрых его изменений. Метод Рунге-Кутты с переменнымшагом разбивает интервал не на равномерные шаги, а более оптимальнымспособом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений — частыми. В результате, длядостижения одинаковой точности требуется меньшее число шагов, чем дляrkfixed.
Метод Булирша-Штера Bulstoer часто оказывается более эффективным для поиска гладких решений.11.3.2. Решение систем ОДУв одной заданной точкеЗачастую при решении дифференциальных уравнений требуется определитьзначения искомых функций не на всем интервале ( t o , t i ) , а только в однойего последней точке. Например, весьма распространены задачи поиска аттракторов динамических систем. Известно, что для широкого класса ОДУодна и та же система при разных (или даже любых, как рассмотренный выше пример осциллятора с затуханием) начальных условиях при t-»°° прихо-278Часть III. Численные методыдит в одну и ту же точку (аттрактор).
Поэтому часто нужно определитьименно эту точку.Такая задача требует меньше ресурсов компьютера, чем решение системыОДУ на всем интервале, поэтому в Mathcad имеются модификации встроенных функций rkadapt и Bulstoer. Они имеют несколько другой наборпараметров и работают быстрее своих аналогов.•rkadapt fyo, to, t i , a c e , D , k , s) — метод Рунге-Кутты с переменным шагом;П bulstoer (yo, to, t i , а с е , D , k , s) — метод Булирша-Штера;• уо — вектор начальных значений в точке to;•t o , t i — начальная и конечная точки расчета;• асе — погрешность вычисления (чем она меньше, тем с лучшей точностью будет найдено решение; рекомендуется выбирать значения погрешности в районе o.ooi);• D — векторная функция, задающая систему ОДУ;• к — максимальное число шагов, на которых численный метод будетнаходить решение;•s — минимально допустимая величина шага.Как легко заметить, вместо числа шагов на интервале интегрирования ОДУ,в этих функциях необходимо задать точность расчета численным методомзначения функций в последней точке.
В этом смысле параметр асе похожна константу TOL, которая влияет на большинство встроенных численцыхалгоритмов Mathcad. Количество шагов и их расположение определяетсячисленным методом автоматически, чтобы обеспечить эту точность. Два последних параметра нужны для того, чтобы пользователь мог искусственноповлиять на разбиение интервала на шаги. Параметр к служит для того,чтобы шагов не было чрезмерно много, причем, нельзя сделать к>юоо. Параметр s — для того чтобы ни один шаг не был слишком малым для появления больших погрешностей при разностной аппроксимации дифференциальных уравнений внутри алгоритма.
Эти параметры следует задавать явно,исходя из свойств конкретной системы ОДУ. Как правило, проведя ряд тестовых расчетов, можно подобрать их оптимальный набор для каждого конкретного случая.Пример использования функции b u l s t o e r для того же примера (11.2—1)приведен в листинге 11.5. В его первых двух строках, как обычно, определяется система уравнений и начальные условия; в следующей строкематрице и присваивается решение, полученное с помощью b u l s t o e r .Структура этой матрицы в точности такая же, как и в случае решениясистемы ОДУ посредством уже рассмотренных нами ранее встроенныхГлава 11. Обыкновенные дифференциальные уравнения279функций (см.
разд. 11.3.1). Однако в данном случае нам интересна толькопоследняя точка интервала. Поскольку сделанное численным методомколичество шагов, т. е. размер матрицы и, заранее неизвестно, то его необходимо предварительно определить. Это сделано в четвертой строкелистинга, присваивающей это число переменной м, в этой же строке оновыведено на экран. В предпоследней строке листинга осуществлен выводрешения системы ОДУ на конце интервала, т. е. в точке t = 50 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующимместом вектора из предыдущей строки).Листинг 11.5. Решение системы двух ОДУУ1D(t , у) :=I-уо - о.
1 • ух J0.1уО :=0Ju := buistoer (уО , 0 , 50 , 10 -5. D , 300 , 0 . 0001М := length (u '*' ) - 150М = 21\7.638 х Ю "2.648 х Ю - 3им,1 = 7.638 х10-3Чтобы попробовать альтернативный численный метод, достаточно в листинге 11.5 заменить ИМЯ фуНКЦИИ buistoer На rkadapt.Внимание!Функции b u i s t o e r и rkadapt (те, что пишутся с маленькой буквы) непредназначены для нахождения решения в промежуточных точках интервала,хотя они и выдают их в матрице-результате. На рис. 11.6 показаны фазовыепортреты рассматриваемой системы ОДУ, полученные с помощью b u i s t o e r(результат листинга 11.5) и с помощью rkadapt (при соответствующей заменетретьей строки листинга 11.5).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
