Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- K = 1e-5; % M
- kd=10; % arbitrary, set to alter convergence speed, numerical stability
- dt = 1e-7; % arbitrary, alter to change convergence speed, numerical stability
- v0 = [0.2; 0.2 ; K]; % M [HA], [A], [X],
- %% for dilute conditions
- %kd=80;
- %dt = 1e-5;
- %v0 = [0.0001; 0.0001 ; K]; % [HA], [A], [X]
- Nstep = 10000;
- t = [0:dt:(Nstep-1)*dt];
- tdilute = t(round(Nstep*[0.1:0.1:0.9]));
- dv_dt = @(kd,K,v) [ -kd*v(1)+kd/K*v(2)*v(3); kd*v(1)-kd/K*v(2)*v(3)] ;
- dv_ = @(v,dt) dv_dt(kd,K,v)*dt;
- dv = @(v) dv_(v,dt);
- v = v0;
- conc = v0;
- for ii=2:Nstep
- vstep = dv(v);
- v = v + [vstep(1); vstep(2); vstep(2)];
- conc(:,ii) = v;
- if any(t(ii) == tdilute)
- v = v/1.2; % <-- change factor to dilute by here
- end
- end
- figure
- subplot(2,1,1)
- plot(t,conc(1,:),'r')
- hold on
- plot(t,conc(2,:),'g--')
- plot(t,conc(3,:),'b--')
- ylabel('[c]')
- legend('[HA]','[A^-]','[H^+]')
- subplot(2,1,2)
- plot(t,-log10(conc(3,:)),'b--')
- ylabel('pH')
- xlabel('time (au)')
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement