2.5.2. Решение систем оду()

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

• rkadapt(y,x1,x2,acc,n,F,k,s) — возвращает матрицу, содержащую таблицу значений решения задачи Коши на интервале от x1 до x2 для системы обыкновенных дифференциальных уравнений, вычисленную методом Рунге-Кутта с переменным шагом и начальными условиями в векторе y (правые части системы записаны в векторе F, n — число шагов, k — максимальное число промежуточных точек решения, и s — минимально допустимый интервал между точками);

• Rkadapt(y,x1,x2,n,F) — возвращает матрицу решений методом Рунге-Кутта с переменным шагом для системы обыкновенных дифференциальных уравнений с начальными условиями в векторе y, правые части которых записаны в символьном векторе F на интервале от x1 до x2 при фиксированном числе шагов n;

• rkfixed(y,x1,x2,n,F) — возвращает матрицу решений методом Рунге-Кутта системы обыкновенных дифференциальных уравнений с начальными условиями в векторе y, правые части которых записаны в символьном векторе F на интервале от x1 до x2 при фиксированном числе шагов n.

Пример 2.25. Решить систему из двух дифференциальных уравнений Ван-дер-Поля, описывающих колебания в нелинейной динамической системе второго порядка, задав начальное и конечное значения x, равными 0 и 20, и используя 100 шагов решения.

Фрагмент документа, представленный на рис. 2.3, решает эту задачу. Заметим, что фактически переменная x здесь есть время, т. е. x=t.

Решение представлено графически в двух видах. Первый вид -это фазовый портрет колебаний, представляющий собой решение y(x) в плоскости (x,y), на которую нанесены все точки (x_n,y_n) для n от 0 до 99 (всего сто точек). Для большей наглядности точки соединены друг с другом отрезками прямых. Фазовый портрет имеет вид закручивающейся спирали и указывает на то, что система генерирует затухающие колебания.

Второй рисунок отражает зависимости y(x) и y’(x) как функции параметра n. Если учесть, что x=t, то можно считать, что на рисунке представлены временные зависимости y(t) и y’(t) в слегка завуалированной форме - шаг времени равен (20-0)/100=20 безразмерным единицам времени. Эти зависимости также подтверждают затухающий характер колебаний.

Рис. 2.3. Решение системы дифференциальных уравнений с применением функции rkfixed

Пример 2.26. Решить самостоятельно представленную выше задачу, используя функцию Rkadapt. Сравните результаты.

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

Если решение системы дифференциальных уравнений имеет вид гладких функций, то вместо описанной ранее функции rkfixed целесообразно применять функцию Bulstoer (y,x1,x2,n,F). Она возвращает матрицу решения системы обыкновенных дифференциальных уравнений, правая часть которых (в виде первых производных неизвестных функций) записана в векторе F(x,y) при заданных в векторе y начальных условиях и при решении на интервале от x1 до x2 для n точек решения, не считая начальной точки. Функция реализует метод Булирша-Штера (Bulirsch-Stoer).