How to describe loop for my ODE equations?

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)

Answer 1

[Bjorn Gustavsson]

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

[Farangis]

Thank you veru much for your great help. :)

It seems that it works properly.

[Bjorn Gustavsson]

Good that it helps!

If the answer solves your problem, then accept it so that other contributors wont waste time looking to answer it.

[Farangis]

Thank you again.

I did it.

[Farangis]

Sorry for bothering you again.

One question: would you please explain this "O0 = [100;rand(99,1)];" to me?

[Farangis]

@Bjorn Gustavsson

[Farangis]

@Bjorn Gustavsson

I want to define initial condition values for O species from j=1 to j=100 in a vector like this O0 = [0.5 0.5 0.5 ... 0.5].

How can I define that vector without writnig 0.5 for 100 times?

[Bjorn Gustavsson]

Ok, my

O0 = [100;rand(99,1)];

was just to make a [100 x 1] array with 100 in the first element and 99 random numbers between 0 and one for the last 99 elements (because I ought to cook up some initial conditions for the example).

If you want 0.5 in all 100 elements you can do something like this:

O0 = 0.5*ones(100,1);

From this question it appears as if you are new to matlab-programming, and as far as I understand the on-ramp material is designed to get people up to speed as effectively as possible, you might consider browsing through that material as you see fit. It helped me increadibly when I realized that matlab has a vast range of handy functions for many tasts from the most trivial to the most specialized, this cuts down programming-time significantly.

[Farangis]

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]

@Bjorn Gustavsson thank you very much for your tip. I will definitely use this on-ramp materials to improve my programming skils.