ODE45 Chemical kinetics

I'm using an ode45 to simulate production of H+ into an environment with time, but H+ just shows up as a horizontal line with time instead of increasing. If I plug in 0 for initial H (Hi) then I don't get a plot at all. Does anyone have any insight as to why H, or y(3) isn't changing over time? I posted my code below

clear all
global p00 
tic
k1 = 8e13; % oxidation of Fe in min^-1*atm^-1*mol^-2*liter^2, Millero 1985 ph>5
k2 = 7% Fly ash MA As describe the experimental leaching kinetics data pseudo second order
%oxidation of Fe2+
Fei = 7; % initial [Fe]
O2i = 7; % intial [O2]
Hi = 1e-5 % initial [H] going to be generated
FAi = 1; % initial [FA]
Asi = 1
p00 = [k1;k2];
y0(1) = Fei;
y0(2) = O2i;
y0(3) = Hi;
y0(4) = FAi;
y0(5) = Asi
dt = 1;
t = (0:dt:1500);
tspan = [0 max(t)];
[T,y] = ode45('RK4odes',tspan,y0);
toc
plot(T,y)
legend('Fe','O2','H','FA','As')
xlabel('time-min')
ylabel('Concentration-M')

function RK4odes1= RK4odes(~,y)
global p00
k1 = p00(1);
k2 = p00(2);
a = 1 %stoichiometric coefficient defining # of moles As reacting with 1 mole of fly ash
kw = 1e-14
RK4odes1(1,1) = -k1.*y(1).*y(2).*(kw./y(3)).^2 %dFe/dt= -k1*Fe*O2*(kw/H)^2
RK4odes1(2,1) = -0.25.*k1.*y(1).*y(2).*(kw./y(3)).^2 %dO/dt = -1/4*k1*(Fe)*(O2)*(kw/H)^2
RK4odes1(3,1) = 2.*k1.*y(1).*y(2).*(kw./y(3)).^2- k2.*(y(4).^0).*y(3)%dH+/dt = 8/4*k1*(Fe)*(O2)*(kw/H)^2-k2*(Flyash-FA)^0*(H)
RK4odes1(4,1) = -k2.*(y(4).^0).*y(3)%dFA/dt = -k(FA)^0*(H+)
RK4odes1(5,1) = (1./a).*k2.*(y(4).^0).*y(3) %dAs/dt=1/a*k2*(FA)*(H+)

Answer 1

[Walter Roberson]

Your equations come out as

     -(Y1*Y2)/(125000000000000*Y3^2)
     -(Y1*Y2)/(500000000000000*Y3^2)
(Y1*Y2)/(62500000000000*Y3^2) - 7*Y3
                               -7*Y3
                                7*Y3

With y0(3) = 0 you get a division by 0 so the third outcome would be infinite the first time, and after that you would have the difference of infinities so you would get nan.

With y0(3) on the order of 1e-5, the (62500000000000*Y3^2) is still a division by a large number, but even so the leading term turns out to be more than -7*Y3. But the balance is close enough that the line is not just horizontal, and is instead cycling rapidly up and down around a value centered on roughly 4.81E-5 . plot(T,y(:,3)) by itself and look more closely at the line especially close to the beginning.

The other lines are not actually straight (though they appear to be on the graphs)... it is just that the slopes are small enough compared to the initial values that the RK methods are comfortable in modeling by a straight line -- when they model as a straight line, the difference is small enough to be within the configured tolerance.

The zoomed in behavior of the third line hints strongly that you should be modeling this as a "stiff" equation.

[Estefania Garcia]

How would I model this as a stiff equation?

[Walter Roberson]

For stiff, change ode45() to ode23s()

You will not see much difference in output, other than that y(:,3) will be much smoother.

tic
k1 = 8e13; % oxidation of Fe in min^-1*atm^-1*mol^-2*liter^2, Millero 1985 ph>5
k2 = 7; % Fly ash MA As describe the experimental leaching kinetics data pseudo second order
%oxidation of Fe2+
Fei = 7; % initial [Fe]
O2i = 7; % intial [O2]
Hi = 1e-5; % initial [H] going to be generated
FAi = 1; % initial [FA]
Asi = 1;
p00 = [k1;k2];
y0(1) = Fei;
y0(2) = O2i;
y0(3) = Hi;
y0(4) = FAi;
y0(5) = Asi;
dt = 1;
t = (0:dt:1500);
tspan = [0 max(t)];
[T,y] = ode23s(@(t,y)RK4odes(t,y,p00),tspan,y0);
toc

Elapsed time is 0.203760 seconds.

plot(T,y)
legend('Fe','O2','H','FA','As')
xlabel('time-min')
ylabel('Concentration-M')

function RK4odes1= RK4odes(~,y,p00)
k1 = p00(1);
k2 = p00(2);
a = 1; %stoichiometric coefficient defining # of moles As reacting with 1 mole of fly ash
kw = 1e-14;
RK4odes1(1,1) = -k1.*y(1).*y(2).*(kw./y(3)).^2; %dFe/dt= -k1*Fe*O2*(kw/H)^2
RK4odes1(2,1) = -0.25.*k1.*y(1).*y(2).*(kw./y(3)).^2; %dO/dt = -1/4*k1*(Fe)*(O2)*(kw/H)^2
RK4odes1(3,1) = 2.*k1.*y(1).*y(2).*(kw./y(3)).^2- k2.*(y(4).^0).*y(3); %dH+/dt = 8/4*k1*(Fe)*(O2)*(kw/H)^2-k2*(Flyash-FA)^0*(H)
RK4odes1(4,1) = -k2.*(y(4).^0).*y(3); %dFA/dt = -k(FA)^0*(H+)
RK4odes1(5,1) = (1./a).*k2.*(y(4).^0).*y(3); %dAs/dt=1/a*k2*(FA)*(H+)
end

[Estefania Garcia]

Thank you!!