Встроенные функции для решения систем ОДУ

В Mathcad 11 имеются три встроенные функции, которые позволяют решать поставленную в форме (2—3) задачу Коши различными численными методами.

• rkfixed(y0, t0, t1, M, D) — метод Рунге-Кутты с фиксированным шагом,
• Rkadapt(y0, t0, t1, M, D) — метод Рунге-Кутты с переменным шагом;
• Buistoer(y0, t0, t1, M, D) — метод Булирша-Штера;
• у0 — вектор начальных значений в точке to размера NXI;
• t0 — начальная точка расчета,
• t1 — конечная точка расчета,
• M — число шагов, на которых численный метод находит решение;
• D — векторная функция размера NXI двух аргументов — скалярного t и векторного у При этом у — искомая векторная функция аргумента t того же размера NXI.

Соблюдайте регистр первой буквы рассматриваемых функций, поскольку это влияет на выбор алгоритма счета, в отличие от многих других встроенных функций Mathcad, например Find=fmd (см разд 11.3.2)

Каждая из приведенных функций выдает решение в виде матрицы размера (M+1)х(N+1). В ее левом столбце находятся значения аргумента t, делящие интервал на равномерные шаги, а в остальных N столбцах — значения искомых функций y₀(t) ,y₁(t), .. ,y_N-1(t), рассчитанные для этих значений аргумента Поскольку всего точек (помимо начальной) м, то строк в матрице решения будет всего M+1

В подавляющем большинстве случаев достаточно использовать первую функцию rkfixed, как это показано в листинге 11.4 на примере решения системы ОДУ модели осциллятора с затуханием (см. разд. 11 2).

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

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

Не забывайте, что векторную функцию D(t,y) следует определять через компоненты вектора у с помощью кнопки нижнего индекса (Subscript) с наборной панели Calculator (Калькулятор) или нажатием клавиши <[>. В третьей строке листинга определено число шагов, на которых рассчитывается решение, а его последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (о, 50).

Как выглядит все решение, показано на рис. 11.3. Размер полученной матрицы будет равен (M+DX(N+I), т.е. Ю1хз. Просмотреть все компоненты матрицы и, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки. Как нетрудно сообразить, на этом рисунке отмечено выделением расчетное значение первого искомого вектора у0 на 12-м шаге u12.1=0.01. Это соответствует, с математической точки зрения, найденному значению У0(6.о)=0.07. Для вывода элементов решения в последней точке интервала используйте выражения типа Uм1=7.523x103.

Рис. 11.3. Матрица решений системы уравнений (листинг 11.4)

Обратите внимание на некоторое разночтение в обозначении индексов вектора начальных условий и матрицы решения В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором столбце — первой компоненты и т. д.

Чтобы построить график решения, надо отложить соответствующие компоненты матрицы решения по координатным осям" значения аргумента и<0> — вдоль оси х, а u<1> и u<2> — вдоль оси У (рис. 11.4). Как известно, решения обыкновенных дифференциальных уравнений часто удобнее изображать не в таком виде, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций. При этом аргумент входит в них лишь параметрически. В рассматриваемом случае двух ОДУ такой график — фазовый портрет системы — является кривой на фазовой плоскости и поэтому особенно нагляден. Он изображен на рис. 11.5 (слева), и можно заметить, что для его построения потребовалось лишь поменять метки осей на u<1> и u<2>, соответственно.

Фазовый портрет типа, изображенного на рис 11.5, имеет одну стационарную точку (аттрактор), на которую "накручивается" решение В теории динамических систем аттрактор такого типа называется фокусом.

Рис. 11.4. График решения системы ОДУ (11.2—1) (листинг 11.4)

Рис. 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 может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо, либо существуют участки медленных и быстрых его изменений. Метод Рунге-Кутты с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений — частыми. В результате, для достижения одинаковой точности требуется меньшее число шагов, чем для rkf ixed. Метод Булирша-Штера Buistoer часто оказывается более эффективным для поиска гладких решений.