Решение систем оду первого порядка

К системам ОДУ первого порядка приходят при описании самых разнообразных задач биологии, физики, химии, химической технологии, экономики и т. д. Для решения систем ОДУ используется метод Рунге-Кутта, последовательно применяемый к каждому дифференциальному уравнению. Однако существуют задачи, которые невозможно решить этим методом, они называются «жесткими». Чаще всего они встречаются при решении задач химической кинетики, в которых скорости отдельных стадий различаются очень значительно, например в миллионы раз. Для их решения используются специальные алгоритмы, так называемые неявные методы. Интересно отметить, что этими методами можно решать и нежесткие, обычные системы ОДУ, однако в жестких методах приходится задавать матрицу частных производных всех функций, что вызывает у некоторых пользователей затруднения и является в некоторых случаях причиной ошибок при формировании матрицы Якоби. Поэтому рекомендуем, сначала попробовать решить задачу простым методом Рунге-Кутта, а в случае неудачи использовать методы решения «жестких» систем.

Решим совместно уравнения (67), описывающие влияние температуры на изменение энтропии и энергии Гиббса в химической реакции. Правые части этой системы являются функцией температуры:

(68)

(69)

Систему дифференциальных уравнений задают в матрице D, куда записывают правые части дифференциальных уравнений (67). В данном случае мы поместили сюда выражения для –_rS⁰_T из уравнения (68) и выражение для _rСр из уравнения (69). Далее в матрицу y0 записываем начальные условия. Решение получаем в виде матрицы U с помощью встроенной функции rkfixed(y0, t0,t1,M,D), где t0 и t1 – начальное и конечное значение аргумента, для нашего случая – температура; М – количество интервалов, на которое разбивается интервал интегрирования, для нашего случая выбрано М=8, что дает вычисление функций при температурах кратных 50 К. Именно для этих температур и можно получить решение.

В программе 32 даны три варианта получения решения. Первый – вывод всей матрицы решения. Видно, что в первом столбце (он считается нулевым) – значения аргумента (у нас – температура), во втором (с номером 1) – значение первой искомой функции (у нас – _rG⁰_T), в третьем (с номером 2) – значение второй искомой функции (у нас – _rS⁰_T). Для построения графиков используются отдельные столбцы матрицы (программа 33).

П рограмма 33