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

В Mathcad имеется несколько встроенных функций, которые позволяют решать задачу Коши различными численными методами. Для "хороших" нежестких систем ОДУ применяются следующие функции:

rkfixed(y0, t0, t1,M,D) — метод Рунге—Кутты с фиксированным шагом;
Rkadapt (y0,t0, t1,M,D) — метод Рунге—Кутты с переменным шагом;
Buistoer(y0,t0,t0,M,D) — метод Булирша—Штера;

у0 — вектор начальных значений в точке t0 размера NXI;

t0 — начальная точка расчета;

t1 — конечная точка расчета;

M — число шагов, на которых численный метод находит решение;

D— векторная функция размера Nx1 двух аргументов — скалярного t и векторного у. При этом у — искомая векторная функция аргумента t того же размера Nx1.

ВНИМАНИЕ!

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

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

В листинге 9.3 приведен пример решения той же самой системы ОДУ осциллятора с затуханием (см. разд. 9.1 и 9.2) при помощи первой из функций rkfixed. Результат расчетов представлен на рис. 9.5 как содержимое матрицы и на рис. 9.6 в виде графика. Точки, в которых получено решение, отмечены на графике рис. 9.6 кружками. Чтобы использовать другой численный алгоритм, достаточно поменять имя функции rkfixed в последней строке листинга на другое (на практике как раз более эффективны функции Rkadapt и Bulstoer).

Листинг 9.3. Решение системы двух ОДУ(модель осциллятора)

Рис. 9.5. Результат, выдаваемый встроенной функцией в качестве решения системы ОДУ (продолжение листинга 9.3}

Первая строка листинга представляет задание параметров модели, вторая — начального условия задачи Коши, а в третьей строке листинга определено число шагов, на которых рассчитывается решение. Самая важная — это предпоследняя строка листинга, в которой, собственно, определяется система ОДУ. Последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (0,40).

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

Матрица, представляющая решение, показана на рис. 9.5. Размер полученной матрицы будет равен (M+i)x(N+l), т. е. 51x3. Просмотреть все компоненты матрицы и, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки. Как нетрудно сообразить, на этом рисунке отмечены выделением расчетные значения искомых векторов на 10-м шаге. Это соответствует, с математической точки зрения, найденным значениям уо(8.0)=о. 307 и yi (8.0) =0.204. Для вычисления элементов решения в последней точке интервала используйте вывод элемента ин, ь Для построения графика решения надо отложить соответствующие компоненты матрицы решения по координатным осям: значения аргумента и<0> — вдоль оси х, а u^<1> и u^<2> — вдоль оси Y (рис. 9.6).

Рис. 9.6. График решения системы ОДУ осциллятора (продолжение листинга 9.3)

ВНИМАНИЕ!

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

Содержание раздела