How to describe loop for my ODE equations?

2 Ansichten (letzte 30 Tage)
Farangis am 8 Feb. 2023
Kommentiert: Farangis am 9 Feb. 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)

Akzeptierte Antwort

Bjorn Gustavsson
Bjorn Gustavsson am 8 Feb. 2023
Bearbeitet: Bjorn Gustavsson am 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 Kommentare
Farangis am 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 am 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

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

Start Hunting!

Translated by