Solving system of simultanous ODE equations with Multiple Initial conditions

2 次查看（过去 30 天）
Meteb Mejbel
Meteb Mejbel 2021-5-5
评论： Alan Stevens 2021-5-7
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.

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2021-5-6
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 个评论
显示 1更早的评论
Meteb Mejbel
Meteb Mejbel 2021-5-6
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
Alan Stevens 2021-5-7
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by