Расчёт концентраций для случая сложного равновесия в Maple
Произвёл попытку решить эту задачу аналитически и численно в Maple:
### A+A<=>B (k1); B+A<=>C (k2); C+A<=>D (k3)
cA:=cA0-x; cB:=k1*cA^2; cC:=k1*k2*cA^3; cD:=k1*k2*k3*cA^4;
### Задаём константы
Conditions:={cA0=1e-40, k1=2.11E+27, k2=7.93E+13, k3=3.48E+10};
### Задаём уравнение материального баланса
MaterialBalance := cA+2*cB+3*cC+4*cD=cA0;
### Выводим его в красивой форме на экран
MaterialBalanceEquation:= collect(expand(lhs(MaterialBalance)-rhs(MaterialBalance)),x)=0;
### Решаем выражение аналитически относительно x
AnlSol:=solve(MaterialBalance,x);
### Решаем частный случай уравнения с подставленными константами
NumSol:=subs(Conditions, AnlSol);
### Вычисляем концентрации
ResultcA:=evalf(subs(Conditions, (cA0 — NumSol)) );
ResultcB:=evalf(subs(Conditions, k1*ResultcA^2));
ResultcC:=evalf(subs(Conditions, k1*k2*ResultcA^3));
ResultcD:=evalf(subs(Conditions, k1*k2*k3*ResultcA^4));
### Расчёт проверочных параметров
SummC:=ResultcA+2*ResultcB+3*ResultcC+4*ResultcD;
Const1:=ResultcB/(ResultcA^2);
Const2:=ResultcC/(ResultcB*ResultcA);
Const3:=ResultcD/(ResultcC*ResultcA);
После всех перестановок получается уравнение:
4*k1*k2*k3*x^4 + (-3*k1*k2-16*k1*k2*k3*cA0)*x^3 +
(2*k1+9*k1*k2*cA0+24*k1*k2*k3*cA0^2)*x^2 +
(-4*k1*cA0-9*k1*k2*cA0^2-1-16*k1*k2*k3*cA0^3)*x +
3*k1*k2*cA0^3+2*k1*cA0^2+4*k1*k2*k3*cA0^4 = 0
Задавая задачу в таком виде получается странное решение.
При нормальных начальных концентрациях c(A0) решение получается нулевым (а нужно посчитать для c(A0)=2.00 и 0.1 ).
При сверхнизких (1.00e-40) всё вроде сходится.