Runge-Kutta for multiple variables

48 次查看（过去 30 天）
Juliana Quintana
Juliana Quintana 2022-4-26
回答： Davide Masiello 2022-4-27
Hello everyone, I am developing a code for resolving a kinetic model of reactions and calculate the concentrations of 4 compounds. I'm implementing the method 4th order Runge-Kutta for multiple variables to calculate each concentration. For this I am coding 4 separate loops for each concentration. I am having trouble with the response vectors; it seems that the course of the reaction isn't correct. I'll appreciate any help solving this issue also an explanation why this is happening will be really helpful to improve my knowledge in Matlab.
clear all;
close all;
clc;
k = [0.75,0.71,0.68];
E0 = 0.5; %Enzima
S0 = 1; %Ssutrato
C0 = 0;%Complejo
P0 = 0; %producto
dt = 10e-4;
x = 0:dt:3 ; %minutos
E = zeros(1,length(x));
S = zeros(1,length(x));
C = zeros(1,length(x));
P = zeros(1,length(x));
fE = @(t,E,S,C)-(k(1).*E.*S)+((k(2)+k(3)).*C);
fS = @(t,E,S,C)-(k(1).*S.*E)+((k(2).*C));
fC = @(t,E,S,C)(k(1).*S.*E)-((k(2)+k(3)).*C);
fP = @(t,C)k(3).*C;
RUNKE-KUTTA 4
%Enzima
E(1) = E0;
for i=1:(length(x)-1)
k_1 = fE(x(i),E(i),S(i),C(i));
k_2 = fE(x(i)+0.5*dt,E(i)+0.5*dt*k_1,S(i),C(i));
k_3 = fE((x(i)+0.5*dt),(E(i)+0.5*dt*k_2),(S(i)),C(i));
k_4 = fE((x(i)+dt),(E(i)+k_3*dt),(S(i)),(C(i)));
E(i+1) = E(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt;
end
%
% %Sustrato
S(1) = S0;
for i=1:(length(x)-1)
k_1 = fS(x(i),E(i),S(i),C(i));
k_2 = fS(x(i)+0.5*dt,E(i),S(i)+0.5*dt*k_1,C(i));
k_3 = fS((x(i)+0.5*dt),E(i),(S(i)+0.5*dt*k_2),C(i));
k_4 = fS((x(i)+dt),E(i),(S(i)+k_3*dt),C(i));
S(i+1) = S(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt;
end
%
% %Complejo
C(1) = C0;
for i=1:(length(x)-1)
k_1 = fC(x(i),C(i),S(i),E(i));
k_2 = fC(x(i)+0.5*dt,E(i),S(i),C(i)+0.5*dt*k_1);
k_3 = fC((x(i)+0.5*dt),E(i),S(i),(C(i)+0.5*dt*k_2));
k_4 = fC((x(i)+dt),E(i),S(i),(C(i)+k_3*dt));
C(i+1) = C(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt;
end
%
% %Producto
P(1) = P0;
for i=1:(length(x)-1)
k_1 = fP(x(i),P(i));
k_2 = fP(x(i)+0.5*dt,P(i)+0.5*dt*k_1);
k_3 = fP((x(i)+0.5*dt),(P(i)+0.5*dt*k_2));
k_4 = fP((x(i)+dt),(P(i)+k_3*dt));
P(i+1) = P(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*dt;
end
figure(1)
plot(x,E,x,S,x,C,x,P)
xlabel('$t (min)$',Interpreter='latex')
ylabel('Concentraciones $(\frac{mmol}{L})$',Interpreter='latex')
legend('Enzima','Sustrato','Compuesto','Producto')
grid on
grid minor

请先登录，再回答此问题。

采纳的回答

Davide Masiello
Davide Masiello 2022-4-27
Ran in:
Hi @Juliana Quintana, the problem is that by defining 4 different function handles and 4 different loops, you have decoupled the ode system.
You can fix this as shown below
clear,clc
k = [0.75,0.71,0.68];
E0 = 0.5; % Enzima
S0 = 1; % Sustrato
C0 = 0; % Complejo
P0 = 0; % Producto
h = 1e-4;
t = 0:h:3 ; % Minutos
y = zeros(4,length(t));
y(:,1) = [E0;S0;C0;P0];
% RUNGE-KUTTA 4
for i = 1:length(t)-1
k1 = odeSystem(t(i),y(:,i),k);
k2 = odeSystem(t(i)+h/2,y(:,i)+h*k1/2,k);
k3 = odeSystem(t(i)+h/2,y(:,i)+h*k2/2,k);
k4 = odeSystem(t(i)+h,y(:,i)+h*k3,k);
y(:,i+1) = y(:,i) + (1/6)*(k1+2*k2+2*k3+k4)*h;
end
E = y(1,:);
S = y(2,:);
C = y(3,:);
P = y(4,:);
figure(1)
plot(t,E,t,S,t,C,t,P)
xlabel('$t (min)$',Interpreter='latex')
ylabel('Concentraciones $(\frac{mmol}{L})$',Interpreter='latex')
legend('Enzima','Sustrato','Compuesto','Producto')
grid on
grid minor
function dydt = odeSystem(t,y,k)
E = y(1);
S = y(2);
C = y(3);
P = y(4);
dydt(1,1) = -(k(1).*E.*S)+((k(2)+k(3)).*C);
dydt(2,1) = -(k(1).*S.*E)+((k(2).*C));
dydt(3,1) = (k(1).*S.*E)-((k(2)+k(3)).*C);
dydt(4,1) = k(3).*C;
end

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by