Plotting output of the ode solver for different initial conditions
3 views (last 30 days)
Show older comments
Hello,
I am trying to solve set of differential equations. I want change one initial condition and plot all derivaties with respect to the varied initial condition at certain point in timespan. I have tried to run loop for it and added deval into code but this doesn't work. Does someone know how this can be done correctly?
For my example I want to calculate all derivatives at t=0.0000001 and I want to plot all derivates with respect to different T values. At t=0.0000001 for all of them.
Thanks.
Quick summary of my problem:
F vector is a fuction of T and time. I want to take derivative of F with respect to time (dF/dt) and calculate value of F at t=0.0000001 and I want to do same calculation for different T values (500:50:1000). Then, I want to plot F vs T graph. Is it possible ?
for T=500:50:1000 %Set of initial condition
[t,F]=ode15s('che515deneme',[0 0.00001],[1*0.25*2*101.3, 0, 0, 0, 0, 0, 0*0.5*2*101.3, 1*0.75*2*101.3, T, 1])
y=zeros(10,1);
y=deval([t,F],0.0000001);
end
subplot(2,1,1)
plot(T,y(:,2),T,y(:,3),T,y(:,4),T,y(:,5),T,y(:,6),T,y(:,10))
xlabel('T(s)');
ylabel('Coverage')
legend('N*','NH*','NH2*','NH3*','H*','*')
subplot(2,1,2)
plot(t,F(:,1),t,F(:,7),t,F(:,8))
xlabel('t(s)');
ylabel('Coverage')
legend('P_N2','P_NH3','P_H2')
The function that I constructed equations
function dFdt = che515deneme(t,F)
dFdt=zeros(10,1);
%F(1)= P_N2
%F(2)= N*
%F(3)= NH*
%F(4)= NH2*
%F(5)= NH3*
%F(6)= H*
%F(7)= P_NH3
%F(8)= P_H2
%F(9)= T
%F(10)= *
A1f=5.6*10^1;
A1r=2.0*10^10;
A2f=6.0*10^13;
A2r=2.8*10^14;
A3f=4.7*10^13;
A3r=1.8*10^13;
A4f=3.3*10^13;
A4r=9.3*10^12;
A5f=5.9*10^13;
A5r=2.1*10^6;
A6f=5.5*10^5;
A6r=2.3*10^13;
Ea1f=33.0*10^3;
Ea1r=137.0*10^3;
Ea2f=86.5*10^3;
Ea2r=41.2*10^3;
Ea3f=60.4*10^3;
Ea3r=8.6*10^3;
Ea4f=17.2*10^3;
Ea4r=64.6*10^3;
Ea5f=83.7*10^3;
Ea5r=0;
Ea6f=0;
Ea6r=89.4*10^3;
R=8.314;
k1f=A1f*exp(-Ea1f/(R*F(9)));
k1r=A1r*exp(-Ea1r/(R*F(9)));
k2f=A2f*exp(-Ea2f/(R*F(9)));
k2r=A2r*exp(-Ea2r/(R*F(9)));
k3f=A3f*exp(-Ea3f/(R*F(9)));
k3r=A3r*exp(-Ea3r/(R*F(9)));
k4f=A4f*exp(-Ea4f/(R*F(9)));
k4r=A4r*exp(-Ea4r/(R*F(9)));
k5f=A5f*exp(-Ea5f/(R*F(9)));
k5r=A5r*exp(-Ea5r/(R*F(9)));
k6f=A6f*exp(-Ea6f/(R*F(9)));
k6r=A6r*exp(-Ea6r/(R*F(9)));
dFdt(1)=k1r*F(2)^2-k1f*F(1)*F(10)^2;
dFdt(2)=k2r*F(3)*F(10)-k2f*F(2)*F(6)-2*k1r*F(2)^2+2*k1f*F(1)*F(10)^2;
dFdt(3)=k2f*F(2)*F(6)+k3r*F(4)*F(10)-k2r*F(3)*F(10)-k3f*F(3)*F(6);
dFdt(4)=k3f*F(3)*F(6)-k3r*F(4)*F(10)-k4f*F(4)*F(6)+k4r*F(5)*F(10);
dFdt(5)=k4f*F(4)*F(6)-k4r*F(5)*F(10)-k5f*F(5)+k5r*F(7)*F(10);
dFdt(6)=2*k6f*F(8)*F(10)^2-2*k6r*F(6)^2-k2f*F(2)*F(6)+k2r*F(3)*F(10)-k3f*F(3)*F(6)+k3r*F(4)*F(10)-k4f*F(4)*F(6)+k4r*F(5)*F(10);
dFdt(7)=k5f*F(5)-k5r*F(7)*F(10);
dFdt(8)=-k6f*F(8)*F(10)^2+k6r*F(6)^2;
dFdt(10)=-2*k1f*F(1)*F(10)^2+2*k1r*F(2)^2+k2f*F(2)*F(6)-k2r*F(3)*F(10)+k3f*F(3)*F(6)-k3r*F(4)*F(10)+k4f*F(4)*F(6)-k4r*F(5)*F(10)+k5f*F(5)-k5r*F(7)*F(10)-2*k6f*F(8)*F(10)^2+2*k6r*F(6)^2;
0 Comments
Accepted Answer
darova
on 25 May 2020
Here is the solution
3 Comments
darova
on 26 May 2020
I understand. Look
T = 500:50:1000; % Set of initial condition
y = zeros(length(T),10);
for i = 1:length(T)
[t,F] = ode15s('che515deneme',[0 0.00001],[1*0.25*2*101.3, 0, 0, 0, 0, 0, 0*0.5*2*101.3, 1*0.75*2*101.3, T(i), 1])
y(i,:) = F(end,:);
end
subplot(2,1,1)
plot(T,y(:,[2:6 10])
xlabel('T(s)');
ylabel('Coverage')
legend('N*','NH*','NH2*','NH3*','H*','*')
subplot(2,1,2)
plot(T,y(:,[1 7 8]))
xlabel('t(s)');
ylabel('Coverage')
legend('P_N2','P_NH3','P_H2')
More Answers (0)
See Also
Categories
Find more on Discrete Data Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!