Решение «жестких» систем оду

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

Несмотря на то, что решение прямых задач химической кинетики – получение зависимости концентрации участников реакции от времени – можно решить обычными способами типа Рунге-Кутта, именно при решении таких задач возникают затруднения. Чтобы этого не случилось все чаще и чаще даже простые системы дифференциальных уравнений решают методами, пригодными для решения жестких систем. Так мы поступим и в нашем курсе.

Для решения «жестких» систем дифференциальных уравнений в МаthCad есть два метода со встроенными функции Stiffb и Stiffr c одинаковыми параметрами (y0,t0,t1,M,D,J), где y0,t0,t1,M,D имеют тот же смысл, что и при использовании метода Рунге-Кутта – матрица начальных условий, начальное и конечное значение аргумента, количество интервалов и вектор правых частей дифференциальных уравнений. Последний параметр J – матрица Якоби, составленная из частных производных правых частей уравнений соответственно по аргументу и по каждой функции. В общем виде для системы:

(70)

(71)

Как получить систему дифференциальных уравнений (70) и и сформировать из нее матрицу Якоби, рассмотрим на примере одного частного случая. Кинетику химической реакции превращения вещества A в С, при этом из А медленно образуется неустойчивое промежуточное вещество В, которое может быстро распадаться как на начальное вещество, так и на продукт реакции (константы скорости отдельных стадий приведены над соответствующими стрелками):

0.1 1000

A  B  C (72)

100

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

• Скорость реакции (изменения концентрации вещества) пропорциональна константе скорости (+, если концентрация увеличивается и – если уменьшается) и концентрации начальных веществ;

• Скорость изменения концентрации вещества складывается (алгебраически) из скоростей реакций, в которых участвует это вещество.

Тогда имеем:

(73)

В соответствии с уравнением (71) для этой системы ОДУ имеем следующую матрицу Якоби:

(74)

Первый столбец в матрице Якоби для решения прямой кинетической задачи всегда состоит из нулей, так как в правой части нет , поэтому частная производная равна нулю.

Обозначив концентрации веществ через [A]=y₀, [B]=y₁, [C]= y₂ нетрудно написать программу для решения этой системы ОДУ и построить графики изменения концентраций участников реакции во времени. Решим эту систему при условии, что в начальный момент времени концентрация вещества А равна 1, а промежуточного и конечного вещества в системе нет.

Программа 34

Из результатов работы программы видно, что концентрация вещества А уменьшается, а С – увеличивается. Концентрация промежуточного вещества очень мала, почти в миллион раз меньше А и С (это как раз и указывает на «жесткость» системы). Концентрация промежуточного вещества сначала увеличивается, а потом падает.