«математические пакеты mathcad и mathematica в решении прикладных химических задач» (стр. 4 из 6)
find(T)=
Полученный ответ и будет корнем уравнения:
Рис. 13 – Решение трансцендентного уравнения в Mathcad
Для решения приведенного трансцендентного уравнения в среде Mathematica можно использовать функцию Solve[]:
Рис. 14 - Решение трансцендентного уравнения в Mathematica
Проверка показывает, что найденное значение температуры найдено верно:
Рис. 15 – Проверка достоверности найденного решения в Mathematica
3.4. Решение системы уравнений
Для решения уравнений в Mathcad применяется блок решения “given – find” применяется и для систем уравнений. В этом случае ответ будет представляться в виде столбца значений и вызываться функцией find(x,y,x…) либо minerr(x,y,x…), где в скобках через запятую перечисляются имена переменных, относительно которых надо решить систему (Рис 16, а).
| x:= 1 y:= 2 | x:= 1 y:= 2 |
 |  |
| а | б |
Рис. 16 – Решение системы уравнений в Mathcad
Можно создать переменные, которым бы мы присвоили значение ответа:
вставить матрицу из n строк и одного столбца (n – количество переменных) и присвоить ей значение minerr(x,y,z…) (Рис 16, б), это делается с помощью раздела меню “Insert > Matrix” или сочетания клавиш “Ctrl+M” [7]
В среде Mathematica воспользуемся функциями: Solve и FindRoot
При использовании функции Solve появляется ошибка, которая говорит о том, что пакет Mathematica не может решить это трансцендентное уравнение алгебраическим методом, т.е. аналитически (This transcendental equation cannot be solved by algebraic methods):
Рис. 17 – Отсутствие аналитического решения для системы уравнений в Mathematica
Следовательно, нужно использовать функцию FindRoot и решать уравнение численно по аналогии с решением в Mathcad:
Рис. 18 – Численное решение системы уравнений в Mathematica
Для проверки полученных корней выполним подстановку в систему и убедимся, что решение верное:
Рис .19 - Проверка достоверности найденных решений системы уравнений в Mathematica
3.5. Дифференциальные уравнения
Численное решение дифференциальных уравнений широко применяется для быстрого и точного решения прямой кинетической задачи: определения зависимости концентраций веществ от времени при заданных начальных концентрациях и скоростях.
Для решения дифференциальных уравнений предусмотрены функции odesolve, rkadapt, rkfixed, bulstoer, а также некоторые другие. Подробную информацию о применимости каждого алгоритма можно найти в разделе Help под названием “differential equation solvers”. Остановимся пока на функции rkfixed. Данная функция интегрирует дифференциальные уравнения и их системы методом Рунге-Кутты четвертого порядка с постоянным (fixed) шагом и может применяться для несложных кинетических схем. К сожалению, эта функция требует специфического ввода условия [[8], [9]].
Пусть имеется простейшая кинетическая схема
Известны значения констант скорости (k1=0,1 с-1, k2=0,05 с-1) и начальных концентраций реагентов (A0=1 М, B0=0 М, C0=0 М). Требуется найти, как будут изменяться концентрации этих веществ во времени в течение 100 секунд.
 |
По закону действующих масс зависимость концентраций реагентов от времени можно записать в виде системы дифференциальных уравнений:
Сначала определим значения констант скорости k1=0,1 и k2=0,05, начальный и конечный моменты времени t0=0 и t1=100, а также число точек N=1000 в интервале [t0, t1], в каждой из которых программа будет искать решение. Затем запишем вектор кинетических уравнений C(t,Y), где Y – массив концентраций, Y0, Y1 и Y2 – концентрации веществ А, В и С соответственно, а также вектор начальных значений C0. Теперь присвоим матрице решений S значение rkfixed (вектор уравнений, начальное значение независимой переменной, конечное значение независимой переменной, число точек, начальные условия) (Рис. 20).
Если теперь вызвать S, то предстанет матрица из 1000 строчек и 4 столбцов. Первый столбец – значения времени, второй, третий и четвертый – соответствующие значения концентраций. Отдельно столбцы можно вызвать, как это показано на Рис.20, с помощью верхнего индекса A:=S<1>. Верхний индекс ставится через “Ctrl+6” или кнопку “Matrix Column” панели “Matrix”.
Рис. 20 – Численное решение системы дифференциальных уравнений с помощью
функции rkfixed в Mathcad
Осталось только построить график (Рис. 21), и прямая кинетическая задача решена:
Рис. 21 – Графическое решение системы дифференциальных уравнений в Mathcad
Функция rkadapt требует такого же синтаксиса, но решает уравнения методом Рунге-Кутты с переменным шагом, анализируя скорость изменения решения. Функция bulstoer реализует метод Булирша-Штера (Bulirsch-Stoer Method) для так называемых жестких систем.
Эти функции удобны в работе и при наличии некоторого опыта затруднений с вводом условия не возникает. Но для тех, кто предпочитает бóльшую наглядность, предусмотрена функция odesolve, работающая в блоке решения.
Решим с ее помощью ту же задачу. Точно так же определите вначале константы скорости и временной интервал. Затем запишем блок решения с тремя дифференциальными уравнениями и граничными условиями (Рис. 22).
Дифференциал
вводится через нажатие клавиш “Shift+/” либо из па-нели инструментов “Calculus”; также можно пользоваться значком производной “ ’ ” (вызывается нажатием “Ctrl+F7”). После того, как записаны все уравнения, вызовите функцию odesolve (вектор переменных, независимая переменная, конечное значение). Полученную зависимость концентрации от времени можно вывести на график либо вызвать в виде столбца значений для переноса в другую программу.
,
Рис. 22 – Решение системы дифференциальных уравнений с помощью функции
odesolve в Mathcad
Выбор метода численного решения также оставляется на усмотрение пользователя. Если кликнуть правой кнопкой мыши на функцию odesolve, появится контекстное меню, в котором можно выбрать метод (с фиксированным и адаптирующимся шагом, а также stiff-метод для жестких систем). Для простых систем результат одинаков, но в случае сложной системы (автоколебания), возможно, придется задуматься о методе интегрирования.
В пакете Mathematica подобные дифференциальные уравнения в численном виде нужно решать с помощью функции NDSolve. В результате получаем интерполяции трех функций для значений t от 0 до 100. С помощью функции Plot[Evaluate[ ]] строим зависимости графически (Рис. 23):
Рис. 23 – Решение системы дифференциальных уравнений в Mathematica
Для проведения физико-комических расчетов часто применяется подход, направленный на упрощение математических моделей, основанный на использовании различных приближений химической кинетики. Одним из наиболее распространенных приближений является метод квазистационарных концентраций. Метод применим для описания многостадийных реакций, на определенных стадиях которых образуются промежуточные вещества (интермедиаты) с высокой реакционной способностью (при кинетическом описании цепных [[10]] и ферментативных [[11]] реакций). Условием применимости метода является малое время жизни интермедиатов по сравнению с временем, за которое состав реакционной системы успевает значительно измениться (условие Христиансена). Принимается, что производные концентраций активных промежуточных частиц по времени равны нулю (теорема Боденштейна).
Таким образом, в системе уравнений, входящих в математическую модель, некоторые дифференциальные уравнения могут быть заменены на алгебраические. Кинетический анализ механизма реакции в этом случае существенно облегчается. Например, механизм термического разложения ацетальдегида CH3CHO в определенных условиях может протекать по цепному механизму, включающему следующие стадии [[12]]: