Clear Filters
Clear Filters

How to describe loop for my ODE equations?

5 views (last 30 days)
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)

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 8 Feb 2023
Edited: Bjorn Gustavsson on 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));
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 Comments
Farangis
Farangis on 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);
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 on 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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!