Как можно повысить точность решения решателей ODE? - PullRequest
0 голосов
/ 07 апреля 2019

Я пытаюсь решить систему дифференциальных уравнений, которые представляют химическую реакцию. 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

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

...