Я пытаюсь решить систему дифференциальных уравнений, которые представляют химическую реакцию. 3 функции (представляющие концентрации) должны иметь значения, которые очень высоки по отношению к изменениям, которые они претерпевают. Таким образом, в данной точке y1 (t1) = 0,1 и dy1 (t2) = 10 ^ -53, что приводит к y1 (t2) = 0,1. Эта ошибка округления приводит к тому, что матрицы, используемые для решения, становятся почти единственными. Я знаю, что обычно эту проблему можно обойти, используя функции digits () и vpa (), однако я не уверен, как заставить решатель ode использовать этот уровень точности.
Я пытался использовать переменную точность на входе, однако решатели ode принимают только входы с двойной точностью. Я пробовал символические решения через dsolve, но нет явного решения. Увеличение относительной терпимости не имеет существенного значения.
n = 10000;
tspan = [0 logspace(-12,-4, n)];
y0 = [0, 0, 0.05, 0.101, 0.1, 0];
options = odeset('RelTol',1e-8,'AbsTol',1e-10);
[T,unper(:,:)] = ode23s(@CombODE, tspan, y0, options);
function dy = CombODE(t, y, flag)
K = [4.45740519775544e-08;194796303866372;4.71e+15;3.8e+16]
Km = [8.46547310572481e+47;1.43069785956689e-49;7.41464500062955e-
40;1.83447707114187e-47]
H = y(1);
O = y(2);
O2 = y(3);
H2 = y(4);
OH = y(5);
H2O = y(6);
R1 = K(1)*H2;
R2 = K(2)*O*O;
R3 = K(3)*O*H;
R4 = K(4)*H*OH;
Rm1 = Km(1)*H*H;
Rm2 = Km(2)*O2;
Rm3 = Km(3)*OH;
Rm4 = Km(4)*H2O;
dy = [2*R1-2*Rm1-R3-R4+Rm3+Rm4;
-2*R2+2*Rm2-R3+Rm3;
R2-Rm2;
-R1+Rm1;
R3-Rm3-R4+Rm4;
R4-Rm4;
];
end
Построенные функции показывают ожидаемую форму, однако в трех функциях, которые начинаются с ненулевых значений, изменение равно нулю. Уравнения являются консервативными, поэтому незначительное увеличение других функций должно проявляться как незначительное уменьшение этих трех функций.