Effacer les filtres
Effacer les filtres

How to describe loop for my ODE equations?

2 vues (au cours des 30 derniers jours)
Farangis
Farangis le 8 Fév 2023
Commenté : Farangis le 9 Fév 2023
I am trying to solve a set of ODE equations with MATLAB ode solver but I do ot know how to define the Loop for my equations. If any one could help me I would be thankful.
These are the ODE equations:
O, I, and S are the molar consentration of components
(N= 100) and (j = 1:100)

Réponse acceptée

Bjorn Gustavsson
Bjorn Gustavsson le 8 Fév 2023
Modifié(e) : Bjorn Gustavsson le 8 Fév 2023
Write a function for these coupled ODEs. Inside that function you can use all the normal programmig tools and tricks you can think of. Something like this should work:
function dOISdt = your_chemistry_ode(t,OIS,kP,kD)
O = OIS(1:end-2);
I = OIS(end-1);
S = OIS(end);
dSdt = kD*I;
dIdt = kP*sum(O(1:end-1)) - kD*I;
dOdt = zeros(size(O(:)))
dOdt(1) = -kP*O(1); % Not exactly clear how you want to treat the first spieces in O
for i1 = 2:numel(O)
dOdt(i1) = kP*(O(i1-1)-kP*O(i1));
end
dOISdt = [dOdt;dIdt;dSdt];
end
Then it should be rather straightforward to integrate this with any of the odeNN-functions:
%% Initial conditions
% I just mock something up, completely without ideas about what chain of
% 100 species decaying towards I and S you model
O0 = [100;rand(99,1)];
I0 = pi;
S0 = 0;
OIS0 = [O0;I0;S0];
kP = 1e6; % reaction-rate, no idea about time-constants either
kD = 1e4;
t_span = linspace(0,1,1e6+1); % same again, steps of ms for 1 s...
[tOut,OISout] = ode45(@(t,IOS) your_chemistry_ode(t,OIS,kP,kD),t_span,OIS0);
HTH
  9 commentaires
Farangis
Farangis le 9 Fév 2023
Hi again,
yes, I am a beginer somehow.
I tried O0 = 0.5*ones(100,1); but I have a problem. It does not take initial conditions for I and S species. Of course I made some changes in code in terms of chemical point of view.
I want the initial conditions to be 0.38 for all of P species and I0 = 0.19 and S0 = 0.17
function dPISdt = Poly_ode(t,PIS,kpeel,kD,DP)
P = PIS(1:198);
I = PIS(199);
S = PIS(200);
dPdt = zeros(size(P(:)));
dPdt(1) = -kpeel*P(1); % the first spieces at i1=1 (P(0)=does not exist)
for i1 = 2:DP-2
dPdt(i1) = kpeel*P(i1-1)-kpeel*P(i1);
end
dIdt = kpeel*sum(P(1:end-1)) - kD*I;
dSdt = kD*I;
dPISdt = [dPdt;dIdt;dSdt];
end
%%%% and here is the ode solver:
DP = 200;
PIS0 = [0.38*ones(1,198) 0.19 0.17];
% reaction-rates [1/min]
kpeel = 0.001;
kD = 0.001;
% time [min]
t_exp = [
48
63
78
93
108
123
138
153
168
181.00
240
300.00
];
[t,PIS] = ode45(@(t,PIS) Poly_ode(t,PIS,kpeel,kD,DP),t_exp,PIS0);
plot (t_exp,PIS(:,1),t_exp,PIS(:,2),t_exp,PIS(:,3));
legend('Polysac','ISA','SHA');
xlabel('time (min)');
ylabel('Concentration [mol/L]');
Farangis
Farangis le 9 Fév 2023
@Bjorn Gustavsson thank you very much for your tip. I will definitely use this on-ramp materials to improve my programming skils.

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by