Solving system of simultanous ODE equations with Multiple Initial conditions

Hi,

I have a system of 3 simultaneous equations (equilibrium reactions). I am trying to solve for the concentrations of three species (CA, CB and CC) which they have multiple initials. So, I have created column vector for each species as shown below

clear 
clc
%input data 
CA0=[5.26211621807567e-14;5.26211621807567e-14;6.87296137179579e-13;
    1.52037833531091e-11;3.99853954838407e-10;1.13330390346014e-08;
    3.31962436509397e-07;9.79450831420185e-06;0.000284543711834628;
    0.00790127173623856;0.198232538062387;3.80176746193761]; %infinite supply of A 
CB0=ones(12,1); %local sites (finite)
CC0=zeros(12,1);
h=0.001;
t=0.1; %longer as possible 
tspan=[0:h:t];
option=odeset('RelTol', 1e-6);
[t, y]= ode15s('kinatic3',tspan, [CA0 CB0 CC0], option);
figure(1)
plot(t,y)
xlabel('time')
ylabel('Concentration') 

and my functions file

function fnc=kinatic3(t, y) 
fnc=zeros(size(y));
CA=y(1);
CB=y(2);
CC=y(3); 
k1=10000;
k2=10;
fnc(1,1)=k2*(CC)^2-k1*CA*CB; 
fnc(2,1)=k2*(CC)^2-k1*CA*CB; 
fnc(3,1)=-k2*(CC)^2+k1*CA*CB; 

So I run the file and get this graph which tells me that nothing happened. I am not sure where the problem is located. I am assuming that the functions should be in a vector form, maybe? I really appreciate your help. Thanks in advance.

Answer 1

[Alan Stevens]

You could try something like:

for i = 1:numel(CA0)
    [t, y]= ode15s(@kinatic3,tspan, [CA0(i) CB0(i) CC0(i)], option);
    figure(1)
    plot(t,y)
    hold on 
end

However, note that you have

fnc(1,1)=k2*(CC)^2-k1*CA*CB; 
fnc(2,1)=k2*(CC)^2-k1*CA*CB; 

These two are identical. Did you mean them to be?

[Meteb Mejbel]

Thanks Alan,

Yes they are the same since I have equilibrium reaction: A+B<-->C. So the rate of A is equal to rate the of B.

[Meteb Mejbel]

It does work. Thank you so much. One more question, once I open y variables at workspace, I just get three columns. How can I see all 36 columns of y and in order such that column 1, 2 and 3 will represent i=1 for CA0(1), CB0(1) and CC0(1), column 4, 5 and 6 will represent i=2 for CA0(2), CB0(2) and CC0(2),....etc?

[Alan Stevens]

Like the following? (Note that the first column is time).

%input data 
CA0=[5.26211621807567e-14;5.26211621807567e-14;6.87296137179579e-13;
    1.52037833531091e-11;3.99853954838407e-10;1.13330390346014e-08;
    3.31962436509397e-07;9.79450831420185e-06;0.000284543711834628;
    0.00790127173623856;0.198232538062387;3.80176746193761]; %infinite supply of A 
CB0=ones(12,1); %local sites (finite)
CC0=zeros(12,1);
h=10^-4;
t=2*10^-3; %longer as possible 
tspan=0:h:t;
option=odeset('RelTol', 1e-6);
Y = tspan';
for i = 1:numel(CA0)
[t, y]= ode15s(@kinatic3,tspan, [CA0(i) CB0(i) CC0(i)], option);
    Y = [Y y];
    figure(1)
    plot(t,y)
    hold on
end
xlabel('time')
ylabel('Concentration') 
disp(Y)

function fnc=kinatic3(~, y) 
fnc=zeros(size(y));
CA=y(1);
CB=y(2);
CC=y(3); 
k1=10000;
k2=10;
fnc(1,1)=k2*(CC)^2-k1*CA*CB; 
fnc(2,1)=k2*(CC)^2-k1*CA*CB; 
fnc(3,1)=-k2*(CC)^2+k1*CA*CB; 
end