Hello Rosalie,
I'm not a chemist so I will take your equations as given although log10 and exp don't seem compatible somehow. Here are a couple of ways to go about this. First a rate function is defined and is at the end of the code, since when running Matlab scripts all functions have to be at the end. The rate function is shown as a function of both time and the amount of Fe. In your case the rate does not explicitly depend on time, only on Fe, but the syntax of the second method requires the time input to be there even if you don't need it. Note that the rate is negative.
I assumed that the listed value for Fe was the starting point at t=0. The first method is simple and just steps through time using the rate calculated from the previous time step, similar to what you had written. The second method uses one of the Matlab differential equation solvers which calculates its own time steps and is quite a bit more accurate. But if you zoom in on the plot you can see that the results are pretty close.
dt = .01;
Fe = zeros(1,1201);
t = zeros(1,1201);
t(1) = 0;
Fe(1) = 0.01*0.26/1000;
for n = 2:1201
t(n) = t(n-1) + dt;
Fe(n) = Fe(n-1) + Ferate(t(n-1),Fe(n-1))*dt;
end
[t2 Fe2] = ode45(@(t,y) Ferate(t,y), [0 12], 0.01*0.26/1000);
figure(1);
plot(t,Fe,t2,Fe2)
legend('method 1','method2')
function rate = Ferate(t,Fe)
Si=1.8/1000;
pH=8;
H=10^(-pH);
k=1.11*10^(-9);
b=1.81;
Ksp=4.68*10^(22);
IAP=(Fe.^(3)*Si^(2))/H^(6);
logIAP=log10(IAP);
logKsp=log10(Ksp);
rate=-k*exp(b*(logIAP-logKsp));
end
Best Answer