Решение систем ОДУ в одной заданной точке

Зачастую при решении дифференциальных уравнений требуется определить значения искомых функций не на всем интервале (to,ti), а только в одной его последней точке. Например, весьма распространены задачи поиска аттракторов динамических систем. Известно, что для широкого класса ОДУ одна и та же система при разных (или даже любых, как рассмотренный выше пример осциллятора с затуханием) начальных условиях приходит в одну и ту же точку (аттрактор). Поэтому часто нужно определить именно эту точку.

Такая задача требует меньше ресурсов компьютера, чем решение системы ОДУ на всем интервале, поэтому в Mathcad имеются модификации встроенных функций rkadapt и Bulstoer. Они имеют несколько другой набор параметров и работают быстрее своих аналогов.

• rkadapt(y0,t0,t1,асе, D, k, s) — метод Рунге-Кутты с переменным шагом;
• buistoer(y0, t0, t1, acc, D, k, s) —метод Булирша-Штера;
• у0 — вектор начальных значений в точке t0;
• t0,t1 — начальная и конечная точки расчета;
• асе — погрешность вычисления (чем она меньше, тем с лучшей точностью будет найдено решение; рекомендуется выбирать значения погрешности в районе 0.001);
• D — векторная функция, задающая систему ОДУ;
• k — максимальное число шагов, на которых численный метод будет находить решение;
• s — минимально допустимая величина шага.

Как легко заметить, вместо числа шагов на интервале интегрирования ОДУ, в этих функциях необходимо задать точность расчета численным методом значения функций в последней точке. В этом смысле параметр асе похож на константу TOL, которая влияет на большинство встроенных численных алгоритмов Mathcad. Количество шагов и их расположение определяется численным методом автоматически, чтобы обеспечить эту точность. Два последних параметра нужны для того, чтобы пользователь мог искусственно повлиять на разбиение интервала на шаги. Параметр k служит для того, чтобы шагов не было чрезмерно много, причем, нельзя сделать k>1000. Параметр s — для того чтобы ни один шаг не был слишком малым для появления больших погрешностей при разностной аппроксимации дифференциальных уравнений внутри алгоритма. Эти параметры следует задавать явно, исходя из свойств конкретной системы ОДУ. Как правило, проведя ряд тестовых расчетов, можно подобрать их оптимальный набор для каждого конкретного случая.

Пример использования функции buistoer для того же примера (11.2—1) приведен в листинге 11.5. В его первых двух строках, как обычно, определяется система уравнений и начальные условия; в следующей строке матрице и присваивается решение, полученное с помощью buistoer. Структура этой матрицы в точности такая же, как и в случае решения системы ОДУ посредством уже рассмотренных нами ранее встроенных

функций (см. разд. 11.3.1). Однако в данном случае нам интересна только последняя точка интервала. Поскольку сделанное численным методом количество шагов, т. е. размер матрицы и, заранее неизвестно, то его необходимо предварительно определить. Это сделано в четвертой строке листинга, присваивающей это число переменной м, в этой же строке оно выведено на экран. В предпоследней строке листинга осуществлен вывод решения системы ОДУ на конце интервала, т. е. в точке t=s0 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующим местом вектора из предыдущей строки).

Листинг 11.5. Решение системы двух ОДУ

Чтобы попробовать альтернативный численный метод, достаточно в листинге 11.5 заменить имя функции buistoer на rkadapt.

Функции buistoer и rkadapt (те, что пишутся с маленькой буквы) не предназначены для нахождения решения в промежуточных точках интервала, хотя они и выдают их в матрице-результате. На рис. 11.6 показаны фазовые портреты рассматриваемой системы ОДУ, полученные с помощью buistoer (результат листинга 11.5) и с помощью rkadapt (при соответствующей замене третьей строки листинга 11.5). Видно, что несмотря на высокую точность (10-5) и верный результат на конце интервала, левый график мало напоминает правильный фазовый портрет (см. рис. 11.5 или правый график на рис. 11.6), начиная быть приемлемым только при предельно допустимом для обсуждаемых функций значении асс=10-16.

В заключение остановимся на влиянии выбора параметра асе на расчеты. Для этого воспользуемся простой программой, представленной на листинre 11.6. В ней из матрицы решения все той же задачи Коши взято лишь полученное значение одной из функций на правой границе интервала. Но зато этот результат оформлен в виде функции пользователя у(е), в качестве аргумента которой выбран параметр асе функции bulstoer.

Рис. 11.6. Фазовый портрет, полученный bulstoer (слева) и rkadapt (справа) (листинг 11.5)

Листинг 11.6. Использование решения ОДУ для определения функции пользователя

Вычисленный вид у(е) показан на рис. 11.7 вместе с аналогичным результатом для функции rkadapt. Как видно, в данном примере численные методы работают несколько по-разному. Метод Рунге-Кутты дает результат тем ближе к истинному, чем меньше выбирается е=асс. Метод Булирша-Штера демонстрирует менее естественную зависимость у (Б): даже при относительно больших е реальная точность остается хорошей (намного лучше метода Рунге-Кутты). Поэтому для экономии времени расчетов (подчеркнем еще раз: для данной конкретной задачи) в функции bulstoer можно выбирать и большие асе.

Чтобы обеспечить заданную точность, алгоритмы, реализованные во встроенных функциях, могут изменять как количество шагов, разбивающих интервал (t0.t1), так и их расположение вдоль интервала. Чтобы выяснить, на сколько шагов разбивался интервал при расчетах у(е)на рис. 11.7 для каждого Е, следует вычислить размер получающейся матрицы. Для этого можно, например, определить функции.

Рис. 11.7. Зависимость расчетного значения одного из уравнений системы ОДУ на конце интервала от параметра асе (листинг 11.6)

Рис. 11.8. Зависимость числа шагов от параметра асе численных методов

Сравнив два результата применения rkadapt для k=30 и k=100, обратите внимание (рис. 11.8), как еще один параметр — максимальное число шагов k, влияет на вид м(е). Заметим, что такие же изменения параметра k на расчет м(е) посредством функции bulstoer влияют слабо.

Таким образом, проводя тестовые расчеты для различных задач и подбирая наилучший набор параметров, можно существенно сэкономить ресурсы компьютера. Конечно, проводить подобный анализ стоит в случаях, когда время расчетов играет важную роль.