Для получения, распределения технологических параметров во времени и в пространстве (в пределах объекта), необходимо произвести СДУ методом, которых дал бы высокую точность решения при минималььных затратах времени на решение, потому что ЭВМ должна работать в режиме реального времени и успевать за ходом технологического процесса. Если время на решение задачи большое, то управляющее воздействие, выработанное на ЭВМ может привести к отрицательным воздействиям. Методов решения существует очень много. В данной работе будет рассмотрен метод решения СДУ методом Рунге-Кутта 4 порядка.

Для удобства работы на ЭВМ, необходимо данную кинетическую схему преобразовать в удобный для работы на компьютере вид. Для этого необходимо кинетическую схему процесса представить в виде уравнений. При рассмотрении кинетической схемы процесса необходимо учитывать коэффициенты скоростей реакций. Но, так как процесс протекает при изотермических условиях, коэффициенты скоростей реакций можно считать за константы скоростей химической реакции. Из приведенной ниже схемы мы можем составить ряд дифференциальных уравнений, учитывающих изотермичность процесса. {2}

Так как коэффициенты K1,K2,K3,K4 являются константами, то можно уравнение записать в следущем виде.{3}

Для преобразования данных дифференциальных уравнений для использования их в расчетах тепловых и кинетических схем методами Рунге-Кутты необходимо подставлять вместо производных значений концентраций, значения концентраций данных в начале процесса. Это обусловлено тем, что метод Рунге-Кутты четвертого порядка, который будет использован для расчета кинетической схемы процесса. Так как этот метод требует сведений только об одной точке и значений функции.

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

Методы Рунге-Кутта обладают следующими свойствами:

1. Эти методы являются одноступенчатыми: чтобы найти у _m+1, нужна

2. Они согласуются с рядом Тейлора вплоть до членов порядка hp, где степень р

3. Они не требуют вычисления производных от f (x, y), а требуют вычисления

Рассмотрим сначала геометрическое построение и выведем некоторые формулы на основе геометрических аналогий. После этого мы подтвердим полученные результаты аналитически.

Предположим, нам известна точка xm, ym на искомой кривой. Тогда мы можем провести прямую линию с тангенсом угла наклона уў _m =f(x _m, y _m ), которая пройдет через точку x _m, y _m. Это построение показано на рис. 1, где кривая представляет собой точное, но конечно неизвестное решение уравнения, а прямая линия L1 построена так, как это только что описано. {4}

Тогда следующей точкой решения можно считать ту, где прямая L ₁ пересечет ординату, проведенную через точку x=x _m+1 =x _m +h.

Уравнение прямой L ₁ выглядит так: y=y _m +yў _m (x-x _m ) так как yў =f(x _m, y _m ) и кроме того, x _m+1 =x _m +h тогда уравнение примет вид

Ошибка при x=x _m+1 показана в виде отрезка е. Очевидно, найденное таким образом приближенное значение согласуется с разложением в ряд Тейлора вплоть до членов порядка h, так что ошибка ограничения равна e _t =Кh2

Заметим, что хотя точка на графике 1 была показана на кривой, в действительности y _m является приближенным значением и не лежит точно на кривой.

Формула 1. 1 описывает метод Эйлера, один из самых старых и широко известных методов численного интегрирования дифференциальных уравнений. Отметим, что метод Эйлера является одним из методов Рунге-Кутта первого порядка.

Рассмотрим исправленный метод Эйлера и модификационный метод Эйлера. В исправленном методе Эйлера мы находим средний тангенс угла наклона касательной для двух точек: x _m, y _m и x _m +h, y _m +hyў _m. Последняя точка есть та самая, которая в методе Эйлера обозначалась x _m+1, y _m+1. Геометрический процесс нахождения точки x _m+1, y _m+1 можно проследить по рис. 2. С помощью метода Эйлера находится точка x _m +h, y _m +hyў _m, лежащая на прямой L _1. В этой точке снова вычисляется тангенс, дает прямую L . Наконец, через точку x _m, y _m мы проводим прямую L, параллельную L . Точка, в которой прямая L пересечется с ординатой, восстановленной из x=x _m+1 =x _m +h, и будет искомой точкой x _m+1, y _m+1.

Тангенс угла наклона прямой L и прямой L равен

Ф(x _m, y _m, h)=Ѕ [f(x _m, y _m )+f(x _m +h, y _m +yў _m h)] 1. 2

где yў _m =f(x _m, y _m ) 1. 3

Уравнение линии L при этом записывается в виде

y=y _m +(x-x _m )Ф(x _m, y _m, h),

так что

y _m+1 =y _m +hФ(x _m, y _m, h). 1. 4

Соотношения 1. 2, 1. 3, 1. 4 описывают исправленный метод Эйлера. {5}

Чтобы выяснить, насколько хорошо этот метод согласуется с разложением в ряд Тейлора, вспомним, что разложение в ряд функции f(x, y) можно записать следующим образом:

f(x, y)=f(x _m, y _m )+(x-x _m )¶ f/¶ x+(y-y _m )¶ f/¶ x+ј 1. 5

где частные производные вычисляются при x=x _m и y=y _m.

Подставляя в формулу 1. 5 x=x _m +h и y=y _m +hyў _m и используя выражение 1. 3 для yў _m, получаем

f(x _m +h, y _m +hyў _m )=f+hf _x +hff _y +O(h2),

где снова функция f и ее производные вычисляются в точке x _m, y _m. Подставляя результат в 1. 2 и производя необходимые преобразования, получаем

Ф(x _m, y _m, h)=f+h/2(f _x +ff _y )+O(h2).

Подставим полученное выражение в 1. 4 и сравним с рядом Тейлора

y _m+1 =y _m +hf+h2/2(f _x +ff _y )+O(h3).

Как видим, исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до членов степени h2, являясь, таким образом, методом Рунге-Кутты второго порядка.

Рассмотрим модификационный метод Эйлера. Рассмотрим рис. 3 где первоначальное построение сделано так же, как и на рис. 2. Но на этот раз мы берем точку, лежащую на пересечении этой прямой и ординатой x=x+h/2. На рисунке эта точка образована через Р, а ее ордината равна y=y _m +(h/2)yў _m. Вычислим тангенс угла наклона касательной в этой точке

Ф(x _m, y _m, h)=f+(x _m +h/2, y _m +h/2*yў _m ), 1. 6

где yў _m =f(x _m, y _m ) 1. 7

Прямая с таким наклоном, проходящая через Р, обозначена через L *. Вслед за тем, мы проводим через точку xm, ym прямую параллельную L *, и обозначаем ее через L _0. Пересечение этой прямой с ординатой x=x _m +h и даст искомую точку x _m+1, y _m+1. Уравнение прямой можно записать в виде y=y _m +(x-x _m )Ф(x _m, y _m, h),

где Ф задается формулой 1. 6. Поэтому

y _m+1 =y _m +hФ(x _m, y _m, h) 1. 8

Соотношения 1. 6, 1. 7, 1. 8 описывают так называемый модификационный метод Эйлера и является еще одним методом Рунге-Кутта второго порядка. Обобщим оба метода. Заметим, что оба метода описываются формулами вида

y _m+1 =y _m +hФ(x _m, y _m, h) 1. 9

и в обоих случаях Ф имеет вид

Ф(x _m, y _m, h)=a ₁ f(x _m, y _m )+a ₂ f(x _m +b ₁ h, y _m +b ₂ hyў _m ), 1. 10

где yў _m =f(x _m, y _m ) 1. 11

В частности, для исправленного метода Эйлера

a ₁ =a ₂ =1/2;

b ₁ =b ₂ =1. {6}

В то время как для модификационного метода Эйлера

a ₁ =0, a ₂ =1,

b ₁ =b ₂ =1/2.

Формулы 1. 9, 1. 10, 1. 11 описывают некоторый метод типа Рунге-Кутты. Посмотрим, какого порядка метод можно рассчитывать получить в лучшем случае и каковы допустимые значения параметров a _1, a _2, b ₁ и b ₂ .

Чтобы получить соответствие ряду Тейлора вплоть до членов степени h, в общем случае достаточно одного параметра. Чтобы получить согласование вплоть до членов степени h2, потребуется еще два параметра, так как необходимо учитывать члены h2f _x и h2ff _y. Так как у нас имеется всего четыре параметра, три из которых потребуются для создания согласования с рядом Тейлора вплоть до членов порядка h2, то самое лучшее, на что здесь можно рассчитывать - это метод второго порядка.

В разложении f(x, y) в ряд 1. 5 в окрестности точки x _m, y _m положим x=x _m +b ₁ h,

Тогда f(x _m +b ₁ h, y _m +b ₂ hf)=f+b ₁ hf _x +b ₂ hff _y +O(h2), где функция и производные в правой части равенства вычислены в точке x _m, y _m.

Тогда 1. 9 можно переписать в виде y _m+1 =y _m +h[a ₁ f+a ₂ f+h(a ₂ b ₁ f _x +a ₂ b ₂ ff _y )]+O(h3).

Сравнив эту формулу с разложением в ряд Тейлора, можно переписать в виде

Если потребовать совпадения членов hf, то a ₁ +a ₂ =1.

Сравнивая члены, содержащие h2f _x, получаем a ₂ b ₁ =1/2.

Сравнивая члены, содержащие h2ff _y, получаем a ₂ b ₂ =1/2.

Так как мы пришли к трем уравнениям для определения четырех неизвестных, то одно из этих неизвестных можно задать произвольно, исключая, может быть, нуль, в зависимости от того, какой параметр взять в качестве произвольного.

Положим, например, a ₂ =w № 0. тогда a ₁ =1-w , b ₁ =b ₂ =1/2w и соотношения 1. 9, 1. 10, 1. 11 сведутся к

y _m+1 =y _m +h[(1-w )f(x _m, y _m )+w f(x _m +h/2w , y _m +h/2w f(x _m, y _m ))]+O(h3) 1. 12

Это наиболее общая форма записи метода Рунге-Кутта второго порядка. При w =1/2 мы получаем исправленный метод Эйлера, при w =1 получаем модификационный метод Эйлера. Для всех w , отличных от нуля, ошибка ограничения равна

e _t =kh3 1. 13

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

y _m+1 =y _m +h/6(R ₁ +2R ₂ +2R ₃ +R ₄ ) 1. 14

где R ₁ =f(x _m, y _m ), 1. 15

R ₂ =f(x _m +h/2, y _m +hR ₁ /2), 1. 16

R ₃ =f(x _m +h/2, y _m +hR ₂ /2), 1. 17

R ₄ =f(x _m +h/2, y _m +hR ₃ /2). 1. 18

Ошибка ограничения для этого метода равна e _t =kh5

так что формулы 1. 14-1. 18 описывают метод четвертого порядка. Заметим, что при использовании этого метода функцию необходимо вычислять четыре раза.

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

• этот метод является одноступенчатым и одношаговым
• требует информацию только об одной точке
• имеет небольшую погрешность
• значение функции рассчитывается при каждом шаге