Solving system of simultanous ODE equations with Multiple Initial conditions

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

Accepted Answer

Alan Stevens
Alan Stevens on 6 May 2021
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?
  3 Comments
Alan Stevens
Alan Stevens on 7 May 2021
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

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!