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);endsubplot(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;
Best Answer