«математические пакеты mathcad и mathematica в решении прикладных химических задач» (стр. 2 из 6)

Одной из основных математических задач, которые приходится решать для таких уравнений, является задача Коши (или начальная задача). Эта задача возникает тогда, когда начальное состояние системы в точке x₀ считается известным (y(x₀)=y₀) и требуется определить ее поведение при x ≠ x₀.

В тех случаях когда невозможно найти аналитическое выражение y(x) (часто так и происходит), применяют численные методы. В основе их построения лежит тот или иной способ замены дифференциального уравнения y'(x) = f(x, y) его дискретным аналогом.

Простейшим и исторически первым численным методом решения задачи Коши является метод Эйлера. В его основе лежит идея графического построения решения дифференциального уравнения. Этот метод дает одновременно и способ нахождения искомой функции в численной (табличной) форме.

Идея метода заключается в том, что промежуток [a, b] разбивается на конечное число n интервалов длиной h точками x₁, x₂,… x_n , и площадь под кривой производной представляется в виде суммы площадей прямоугольников

Таким образом, на малом промежутке h изменения независимой переменной x_i <x<x_i+h = x_i+1 вместо интегральной кривой дифференциального уравнения y(x) берется касательная к ней

В результате неизвестная интегральная кривая заменяется приближенной к ней ломаной линией (ломаной Эйлера), для которой угловой коэффициент n-го звена равен f(x_n,y_n).

Метод Эйлера обладает невысокой точностью, имея лишь первый порядок аппроксимации: с изменением h ошибка вычислений также меняется пропорционально h. К тому же погрешность каждого нового шага, вообще говоря, систематически возрастает. Наиболее приемлемым для практики методом оценки точности является в данном случае метод двойного счета - с шагом h и шагом h/2. Совпадение полученных двумя способами результатов дает естественные основания считать их верными.

Гораздо более точным является метод Рунге-Кутты. Формулы для наиболее популярного алгоритма четвертого порядка имеют вид:

где коэффициенты k вычисляются как:

Для повышения точности расчета иногда применяется метод Рунге-Кутты с переменным шагом. В этом случае интервал разбивается на участки разной длины: там, где решение меняется слабо, шаги выбираются редкими, а в областях его сильных изменений - частыми. Это очень просто осуществить, так как алгоритм Рунге-Кутты является одношаговым и подразумевает простой пересчет при любом значении шага h_i искомого y(x_i+h_i) через y(x_i). В результате применения адаптированного алгоритма для достижения одинаковой точности может потребоваться существенно меньшее число шагов, чем для стандартного метода Рунге-Кутты с фиксированным шагом.

Несмотря на кажущуюся универсальность, метод Рунге-Кутты «не работает» в случае так называемых жестких дифференциальных уравнений и соответствующих систем. Будьте осторожны в выборе метода, когда решаете, например, уравнение вида y'(x)=-25y + cos(x) + 25sin(x), y(0) = 1, или анализируете кинетическую схему с сильно различающимися константами скорости, типичную для автокатализа. Если вместо красивых кинетических кривых метод Рунге-Кутты выдает непонятные осцилляции, знайте: Вы столкнулись со случаем жесткой системы. Для решения требуется применение совсем иных алгоритмов, например неявных методов типа метода Булирша-Штера.

Основной идеей метода является вычисление состояния системы в точке x+h как результата двух шагов длины h/2, четырех шагов длины h/4, восьми шагов длины h/8 и так далее с последующей экстраполяцией результатов. Метод строит рациональную интерполирующую функцию, которая в точке h/2 проходит через состояние системы после двух таких шагов, в точке h/4 проходит через состояние системы после четырех таких шагов, и т. д., а затем вычисляет значение этой функции в точке h = 0, проводя экстраполяцию. Таким образом, проводится один шаг метода, после чего принимается решение, следует ли изменять шаг, а если да, то в какую сторону. При этом используется оценка погрешности, которую мы получаем в качестве дополнительного результата при рациональной экстраполяции [[3]].

2.2. Гармонический анализ и ряд Фурье

Теорема Фурье гласит, что любую зависящую от времени функцию можно представить в виде суммы синусов и косинусов с соответствующими весами (амплитудами). Таким образом, из общего «хора» сигналов преобразование Фурье выделяет отдельные «голоса», определяет их частоту и амплитуду. Математически можно сказать, что любая периодическая зависимость при определенных (и справедливых для любых реальных сигналов) допущениях может быть представлена бесконечным рядом Фурье. В тригонометрической форме он имеет вид:

где a_k и b_k - косинусные и синусные коэффициенты Фурье. Они вычисляются по формулам:

В этих выражениях k-номер гармоники, f₁ – частота первой гармоники, T= 1/f₁ – период колебаний.

Непериодические сигналы также могут иметь ряды Фурье. Но вместо дискретного спектра (из отдельных гармоник с номерами k=1,2,3,...) они имеют сплошной спектр [[4]].

Применение преобразования Фурье произвело настоящий переворот в спектроскопии. Стало возможным регистрировать и обрабатывать ЯМР-, ИК- и УФ-спектры за несколько секунд или даже до-лей секунд, обеспечивая быстрый перевод первичного сигнала к привычному виду (Рис. 1). Конструкция приборов упростилась, а круг возможностей спектральных инструментов зачастую определяется вычислительной мощностью компьютера [[5]].

Рис. 1 - Сигнал спада свободной индукции (а) и полученный из него спектр ЯМР (б)

Смысл преобразования Фурье заключается в представлении функции в виде суммы периодических функций (синусоид). Основное применение преобразования Фурье в химии - перевод функции от времени в функцию от частоты:

Для численных расчетов был разработан быстрый и эффективный метод быстрого преобразования Фурье (БФП, в англоязычной литературе FFT - Fast Fourier Transform,), включенный во все серьезные математические программы. Смысл БФП заключается в разбиении массива точек на два подмассива и проведении преобразования Фурье над каждым в отдельности. Единственным ограничением БПФ является то, что длина набора данных должна быть равна целой степени двойки (т. е. можно обработать 512 точек, но не 511 или 513) [3, 5].

3. РЕШЕНИЕ НЕКОТОРЫХ ХИМИЧЕСКИХ ЗАДАЧИС ПОМОЩЬЮ МАТЕМАТИЧЕСКИХ ПАКЕТОВ MATHCAD И MATHEMATICA

Система Mathcad уже довольно длительное время является бесспорным лидером среди математического ПО с WYSIWYG (what you see is what you get) интерфейсом, максимально приближающим внешний вид документов к традиционным расчетам «на бумаге». Язык ввода формул прост и понятен интуитивно. К одному из главных отличий и достоинств пакета Mathematica следует отнести возможности получать решения сложных систем уравнений и функций в аналитическом виде, что не под силу Mathcad. Однако, работа в Mathematica требует не только высокого уровня владения пользователем лексики приложения, но и внимательного выполнения поставленной задачи.

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

3.1. Простые расчеты в Mathcad и Mathematica

Следует отметить, что изначально система Mathcad позиционировалась разработчиками как суперкалькулятор. Для примера использования данной системы в простых расчетах, попробуем посчитать какое-нибудь выражение, например, оценить константу равновесия для реакции паровой конверсии СО (ΔH = –41,17 Дж/моль, ΔS = –42,09 Дж/моль*К) при температуре 600 K. Энтальпию и энтропию сорбции в первом приближении можно считать не зависящими от температуры.

Вводим следующее выражение:

exp(-(-41170+42.09*600)/600*8.314)=

На экране появится:

Рис. 2 - Вычисления в Mathcad

Как видно, форма записи максимально приближена к вычислениям «на бумаге».

Численный ответ выдается в виде десятичной дроби, с тремя знаками после запятой. Эти параметры можно изменить в меню “Format > Result”. Представить ответ в виде натуральной дроби позволяет функция Fraction.

Несмотря на то, что ядро Mathcad предназначен для численного решения, он позволяет также производить и несложные символьные расчеты. Например, брать неопределенные интегралы типа интеграла зависимости теплоемкости от температуры a+bT+c/T². Для этого следует ввести знак неопределенного интеграла (Ctrl+I либо с панели “Calculus”), записать уравнение и поставить после него вместо обычного знака равенства значок символьного решения “→” (Ctrl+. либо с панели “Evaluation” [[6], [7]]: