Встроенные функции для решения систем ОДУ
Обратите внимание на некоторое разночтение в обозначении индексов вектора начальных условий и матрицы решения. В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором столбце – первой компоненты и т. д.
Чтобы построить график решения, надо отложить соответствующие компоненты матрицы решения по координатным осям" значения аргумента u<0> – вдоль оси х, а u<1> и u<2> – вдоль оси У (рис. 11.4). Как известно, решения обыкновенных дифференциальных уравнений часто удобнее изображать не в таком виде, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций. При этом аргумент входит в них лишь параметрически. В рассматриваемом случае двух ОДУ такой график – фазовый портрет системы – является кривой на фазовой плоскости и поэтому особенно нагляден. Он изображен на рис. 11.5 (слева), и можно заметить, что для его построения потребовалось лишь поменять метки осей на u<1> и u<2>, соответственно.
Рис. 11.4. График решения системы ОДУ (11.2-1) (листинг 11.4)
Фазовый портрет типа, изображенного на рис 11.5, имеет одну стационарную точку (аттрактор), на которую "накручивается" решение В теории динамических систем аттрактор такого типа называется фокусом.
Рис. 11.5. Фазовый портрет решения системы ОДУ при M=100 (слева) и М=200 (справа) (листинг 11.4)
В общем случае, если система состоит из N ОДУ, то фазовое пространство является N-мерным. При N>3 наглядность теряется, и для визуализации фазового портрета приходится строить его различные проекции.
На том же рис. 11.5, справа, показан для сравнения результат расчета фазового портрета с большим числом шагов. Видно, что в этом случае обеспечивается лучшая точность и, в результате, решение получается более гладким. Конечно, при этом увеличивается и время расчетов.
При решении систем ОДУ многие проблемы могут быть устранены при помощи простой попытки увеличить число шагов численного метода. В частности, сделайте так при неожиданном возникновении ошибки "Found a number with a magnitude greater than 10^307" (Найдено число, превышающее значение 10307). Данная ошибка может означать не то, что решение в действительности расходится, а просто недостаток шагов для корректной работы численного алгоритма.
В заключение следует сказать несколько слов об особенностях различных численных методов. Все они основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации, получаются алгоритмы различной точности и быстродействия. В Mathcad использован наиболее популярный алгоритм Рунге-Кутты четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо rkfixed другие функции, благо сделать это очень просто, благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе.
Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо, либо существуют участки медленных и быстрых его изменений. Метод Рунге-Кутты с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений – частыми. В результате, для достижения одинаковой точности требуется меньшее число шагов, чем для rkfixed. Метод Булирша-Штера (Buistoer) часто оказывается более эффективным для поиска гладких решений.