How to describe loop for my ODE equations?

1 visualizzazione (ultimi 30 giorni)
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)

Risposta accettata

Bjorn Gustavsson
Bjorn Gustavsson il 8 Feb 2023
Modificato: Bjorn Gustavsson il 8 Feb 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));
dOISdt = [dOdt;dIdt;dSdt];
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);
  9 Commenti
Farangis il 9 Feb 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);
dIdt = kpeel*sum(P(1:end-1)) - kD*I;
dSdt = kD*I;
dPISdt = [dPdt;dIdt;dSdt];
%%%% 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 = [
[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));
xlabel('time (min)');
ylabel('Concentration [mol/L]');
Farangis il 9 Feb 2023
@Bjorn Gustavsson thank you very much for your tip. I will definitely use this on-ramp materials to improve my programming skils.

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by