Математические задачи в пакете MathCAD 12
Глава 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)
ВНИМАНИЕ!
Обратите внимание на некоторое разночтение в обозначении индексов вектора начальных условий и матрицы решения. В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором столбце — первой компоненты и т. д.
