Используется для повышения точности получаемого разностного решения задачи Коши (5.1)-(5.2). Расчетные формулы разностной схемы имеют несколько иной вид. Разностное уравнение:
(6.14)
g1,k = f(xk,yk
g2,k = f(xk + h/2,yk + h/2*g1,k) (6.15)
g3,k = f(xk + h/2,yk + h/2*g2,k)
g4,k = f(xk + h,yk + h*g3,k)
Формулы (5.15) обеспечивают точность четвертого порядка относительно h как в аппроксимации, так и в погрешности решения.
Формулу (5.14) можно переписать следующим образом:
yk+1 = yk + h/6(g1,k + 2g2,k + 2g3,k + g4,k) (6.16)
Пример 6.3.
Определить решение дифференциального уравнения
методом Рунге-Кутты 4-ого порядка точности на отрезке [0,1] с шагом h=0,5 и начальным условием: u(0)= - 0,5. Сравнить его с частным аналитическим решением уравнения. Найти погрешность решения z.
Решение.
Его общее решение имеет вид: .
Определим решение методом Эйлера с шагом 0,5:
y(0)=u(0)= - 0,5
g1,0 = f(x0,y0)= 0/( - (-0,5))=0
g2,0 = f(x0 + h/2,y0 + h/2*g1,0) = 0,25/( -(-0.5+0.25*0))= 0,236068
g3,0 = f(x0 + h/2,y0 + h/2*g2,0)= 0,25/( - (-0.5+0.25*0.236068))= 0,263741
g4,0 = f(x0 + h,y0 + h*g3,0)= 0,5/( - (-0.5+0.25*0.263741))= 0,505545
y(0,5) = y(0)+h( g1,0 + 2g2,0 +2 g3,0 + g4,0 )/6= -0,37457
g1,1 = 0,500344
g2,1 = 0,72123
g3,1 = 0,773984
g4,1 = 1,012499
y(1,0) = y(0,5)+h( g1,1 + 2g2,1 +2 g3,1 + g4,1 )/6= 0,000703
Расчетная таблица имеет вид:
x | g1 | g2 | g3 | g4 | y_R-K | u_analit | z | Psi |
0 | 0 | 0,236 | 0,264 | 0,506 | -0,5 | -0,5 | 0 | 0,25 |
0,5 | 0,500 | 0,721 | 0,774 | 1,012 | -0,375 | -0,375 | 0,0004 | 0,250 |
1 | 0,001 | 0 | 0,0007 |
Аналитическое решение дает значения: u(0)=-0,5; u(0,5)= -0,375; u(1,0)= 0
Z = max|zk| = max{0,0.00043, 0.000703}=0.000703
y = max|yk| = max{0.25; 0.249656}=0.25
Соответствие метода Рунге-Кутта и аналитического решения показано на рисунке.
По нему видно, что аналитическое решение очень хорошо согласуется с полученными расчетными значениями, то есть расчетная кривая близка к аналитическому решению и, поэтому, иногда за аналитическое решение принимается кривая, полученная методом Рунге-Кутты или, вообще говоря, методом высшего порядка.
Пример 6.4.
Дана кинетическая схема процесса
A + B (k1)àC
C (k2)àА
Начальные условия: A0=0,03; B0=0,04; C0=0,01, k1=4, k2=3
а) выписать систему дифференциальных уравнений, каждое из которых представляет собой скорость изменения концентрации отдельного вещества в течение процесса;
б) найти решение системы методом Эйлера 1-ого порядка точности на отрезке по времени [0,3] с шагом 0,5;
в) найти решение системы методом Рунге-Кутта 4-ого порядка точности на отрезке по времени [0,3] с шагом 0,5;
г) сравнить кривые зависимости концентрации вещества А от времени, полученные двумя методами.
Решение.
а) По этой кинетической схеме можно составить систему дифференциальных уравнений:
С начальными данными: A0=0,03; B0=0,04; C0=0,01.
б) Все расчеты, связанные с методом Эйлера удобно вести в пакете Excel. Для этого необходимо составить расчетную таблицу. Оформим формулы, соответствующие расчету значений всех функций при времени t=0,5:
B3 | = B2+$A$3*(-4*$B2*$C2+3*$D2)
C3 | = C2+$A$3*(-4*$B2*$C2)
D3 | = D2+$A$3*(4*$B2*$C2-3*$D2)
Однако, такой расчет дает значение в ячейке D3 отрицательное число, что противоречит условию неотрицательности концентрации, поэтому формулу в ячейке D3 необходимо заменить на следующую:
D3 | = ЕСЛИ(D2+$A$3*(4*$B2*$C2-3*$D2)<0;0;D2+$A$3* (4*$B2* $C2 -3*$D2))
Это изменение позволит исправить ситуацию с отрицательными значениями концентраций веществ. Вообще говоря, условие неотрицательности концентраций могут быть использованы для всех веществ: А, В и С.
Тогда таблица Excel будет иметь вид:
Если провести анализ полученных значений и, собственно, системы дифференциальных уравнений, приведенной в пункте а), то видно, что с увеличением времени процесса, концентрация вещества В постепенно уменьшается, об этом же и свидетельствует дифференциальное уравнение dB/dt, оно всюду отрицательно, т.е. отрицательна скорость изменения концентрации вещества, следовательно, оно уменьшается.
Вещества А и С ведут себя с точностью до наоборот, что также видно по их дифференциальным уравнениям. Для этих веществ можно выписать балансовое уравнение:
A(t)+C(t) = A(0) + C(0) = const
Если вновь посмотреть на систему дифференциальных уравнений, то можно заметить, что
dA/dt + dC/dt = 0
то есть суммарная концентрация веществ А и С постоянна и неизменна во времени. Фактически это означает, что если концентрация вещества А увеличивается, то падает концентрация вещества С, и наоборот.
Приведем графики изменения концентраций всех веществ в одной координатной плоскости:
Рисунок свидетельствует обо всех выводах, полученных путем анализа системы дифференциальных уравнений, описывающих заданную кинетическую схему процесса.
Осью симметрии для кривых, характеризующих концентрации веществ А и С, является их средняя суммарная концентрация, равная:
s = (A(0) + C(0))/2 = 0,02
в) Составим расчетные таблицы, для этого внесем в ячейки B3, C3, D3, E3,F3 следующие формулы:
B3 | = -4*$F2*$F11 + 3*$F20
С3|=-4*($F2+0,5/2*B3)*($F11+0,5/2*B12)+3*($F20+0,5/2*B21)
D3|=-4*($F2+0,5/2*C3)*($F11+0,5/2*C12)+3*($F20+0,5/2*C21)
E3|=-4*($F2+0,5*D3)*($F11+0,5*D12)+3*($F20+0,5*D21)
F3|=F2+0,5/6*(B3+2*C3+2*D3+E3)
Аналогично будут выглядеть формулы в ячейках, используемых для расчета коэффициентов вещества В – это ячейки B12, C12, D12, E12 и F12:
B12 | = -4*$F2*$F11
С12|=-4*($F2+0,5/2*B3)*($F11+0,5/2*B12)
D12|=-4*($F2+0,5/2*C3)*($F11+0,5/2*C12)
E12|=-4*($F2+0,5*D3)*($F11+0,5*D12)
F12|=F11+0,5/6*(B12+2*C12+2*D12+E12)
Формулы, находящиеся в ячейках B21, C21, D21, E21, F21:
B21 | = 4*$F2*$F11 - 3*$F20
С21|=4*($F2+0,5/2*B3)*($F11+0,5/2*B12)-3*($F20+0,5/2*B21)
D21|=4*($F2+0,5/2*C3)*($F11+0,5/2*C12)-3*($F20+0,5/2*C21)
E21|=4*($F2+0,5*D3)*($F11+0,5*D12)-3*($F20+0,5*D21)
F21|=F20+0,5/6*(B21+2*C21+2*D21+E21)
Таблица в Excel будет иметь вид, показанный на нижеследующем рисунке. В ней для расчета значений каждого вещества используются свои коэффициенты g1, g2, g3, g4.
По расчетным таблицам можно построить графики изменения концентраций веществ А, В и С в одной координатной плоскости:
Полученный рисунок также свидетельствует об адекватности всех полученных выводов из анализа системы дифференциальных уравнений. Явно прослеживается симметрия между графиками, характеризующими изменение концентрации веществ А и С, причем осью симметрии является число
s = (A(0) + C(0))/2 = 0,02
также как при использовании метода Эйлера.
Ниспадающий вид кривой по изменению концентрации вещества В аналогичен полученному выше методом Эйлера.
Можно сделать вывод, что полученные численными методами зависимости изменения концентраций участвующих в процессе веществ внешне схожи друг с другом, но разнятся значениями.
Так, можно отметить, что наблюдавшаяся отрицательность концентрации вещества С при использовании метода Эйлера устранена с использованием метода Рунге-Кутты. Это является следствием высокой точности последнего метода и можно утверждать, что именно такое поведение концентрации вещества С допустимо в данном случае.
г) Сравнение методов.
Прежде чем перейти непосредственно к построению зависимости изменения концентрации вещества А, полученной двумя численными методами, построим в одной координатной плоскости решения примеров 6.1 и 6.3 совместно.
Как видно по рисунку, метод Эйлера не точно аппроксимирует аналитическое решение задачи (поскольку он первого порядка точности), а метод Рунге-Кутты очень хорошо согласуется с аналитическим решением (его точность четвертого порядка).
Этот результат является важным доказательством того, что при исследовании сложных химических процессов получение необходимых зависимостей лучше вести методом Рунге-Кутты.
Дата: 2019-03-05, просмотров: 273.