龙格库塔4阶求解酶动力学方程

RK4求解酶动力学方程，由两个部分构成：主函数RK4用于进行迭代求解，子函数f用于存储微分方程。

function RK4

clear;clc;
% define initial values
E0=1;  
k3=150/60;   

% define step width and initial values of C,S
h=0.001;        % step width is 0.001s
t=0:h:30;       % range is 0-30s
n=length(t);   % quantity of t
Y(1,1)=0;     % initial values of C put in matrix Y
Y(2,1)=10;      % initial values of S put in matrix Y

% define RK4 to find out the solution of S and C
for k=1:n-1  
    z1=f(t(k),Y(1:2,k));  % Y(1:2,k) means take the 1st and 2nd row of kth line
    z2=f(t(k)+h/2,Y(1:2,k)+z1*h/2);
    z3=f(t(k)+h/2,Y(1:2,k)+z2*h/2);
    z4=f(t(k)+h,Y(1:2,k)+z3*h);
    Y(1:2,k+1)=Y(1:2,k)+h*(z1+2*z2+2*z3+z4)/6; % iteration of Y 
                                               % new C/S will be added into matrix Y
                                               % C to the 1st row, S to the 2nd row
end
C=Y(1,:);  % Y's 1st row
S=Y(2,:);  % Y's 2nd row

%solution to E and P
E=E0-C;

for i=1:1:n  % initial value: step width: final value
    P(i)=k3*sum(C(1:i))*h;  % P is the integration of C, C is discrete
                            % P(0) = 0
                            % so Pi is the sumption of C1-Ci * step wid
end

Vp=k3.*C; % Vp is the change rate of P
          % n2.means conduct every element in matrix

% plot
figure(1);  % 1st picture
plot(t,E,t,S,t,C,t,P,'LineWidth',3);  % x:t y:E,S,C,P
legend('E','S','C','P');  
title('Variation of Components Concentration');
xlabel('time/s');
ylabel('concentration/uM');

figure(2);  % 2nd picture
plot(t,Vp,t,S,'LineWidth',3);  % x:t y:Vp,S
legend('Vp','S');
title('Variation of Vp and S');
xlabel('time/s');

figure(3);  % 3rd picture
plot(S,Vp,'LineWidth',3);  % x:s y:Vp
title('Vp changed with the concentration of S');
xlabel('Sconcentration/uM');
ylabel('uM/s');

end

function F=f(t,Y)  
% input t and matrix Y
% define variables
k1=100/60;
k2=600/60;
k3=150/60;
E0=1;
theta=k2+k3;
lamda=k1*E0;

% define the equation we want to solve
C=Y(1,1);  % Y is a 1*2 matrix
S=Y(2,1); 
f1=lamda*S-(theta+k1*S)*C;  % equation of dC/dt
f2=-lamda*S+(k1*S+k2)*C;  % equation of dS/dt
F=[f1;f2];  % output of F is a matrix contain f1(C) and f2(S)
end

最后用RK4方法求解出的效果如图：