Hi everyone,
I am trying to solve a system of 16 differential equations,until the 14 differential equations I am able to get the results perfectly but adding 15 and 16 equation i am getting the Warning: Failure at t=8.801889e-07. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.694066e-21) at time t.Please help me to solve this
function dy=adia_real(t,y)
k=deadkineticestimation(y(15));
%deadkineticestimation function given below
dy=zeros(16,1);
%equation 15
Hra=50.100;%%j/mol
Hrp=-23.300;
ra=((k(1)*y(1)*y(9))-(k(2)*y(4)*y(6)));
rp=((k(3)*y(5)*y(6))-(k(4)*y(6)));
mtot=11.704;
mtot1=11704;
cp=2.06;
Mmon=144.14;
Mcat=405.1;
Mcocat=186.34;
Macid=144.22;
%equation 16
density_cat=1200;
density_cocat=820;
density_acid=908.8;
density_pla=k(17);
rmon=-rp;
rcat=-ra;
rcocat=((k(2)*y(4)*y(6))-(k(1)*y(1)*y(9)));
racid=ra;
a=rmon/k(17);
b=(y(5)*y(16)*Mmon)/((k(17))^2);
diff_mon=((-0.8462)/(1.110865+(7.391*10^-4*y(15)))^2);
diff_pla=((-0.8462)/(1.110865+(7.391*10^-4*y(15)))^2);
c=rcat/density_cat;
d=rcocat/density_cocat;
e=racid/density_acid;
f=(1-(1/density_pla));
g=(mtot-(((y(5)*Mmon)+(y(1)*Mcat)+(y(2)*Mcocat)+(y(4)*Macid))*y(16))/(density_pla)^2);
if(y(6)==0.0 || y(7)==0.0 || y(9)==0.0 || y(10)==0.0 || y(12)==0.0 || y(13)==0.0)
lam3=0.0;
mu3=0.0;
gam3=0.0;
else
lam3=(y(8)*((2*y(8)*y(6))-(y(7)*y(7))))./(y(7)*y(6));
mu3=(y(11)*((2*y(11)*y(9))-(y(10)*y(10))))./(y(10)*y(9));
gam3=(y(14)*((2*y(14)*y(12))-(y(13)*y(13))))./(y(13)*y(12));
end
dy(1)=(-1*k(1)*y(1)*y(2))+(k(2)*y(3)*y(4))-(k(3)*y(1)*y(9))+(k(4)*y(6)*y(4))+(rcat-(y(1)/y(16))*dy(16));
dy(2)=(-1*k(1)*y(1)*y(2))+(k(2)*y(3)*y(4))-(k(9)*y(2)*y(6))+(k(10)*y(3)*y(9));
dy(3)=(-1*k(5)*y(5)*y(3))+(k(1)*y(1)*y(2))-(k(2)*y(3)*y(4))+(k(9)*y(2)*y(6))-(k(10)*y(9)*y(3)); % +(k(4)*y(6));
dy(4)=(k(1)*y(1)*y(2))-(k(2)*y(3)*y(4))+(k(3)*y(1)*y(9))-(k(4)*y(6)*y(4)+(racid-(y(4)/y(16))*dy(16)));
dy(5)=(-1*k(5)*y(5)*y(3))-(k(7)*y(5)*y(6))+(k(8)*y(6)+(rmon-(y(5)/y(16))*dy(16)));
dy(6)=(k(3)*y(9)*y(1))-(k(4)*y(6)*y(4))+(k(5)*y(5)*y(3))-(k(9)*y(6)*y(2))+(k(10)*y(9)*y(3)+(rmon-(y(6)/y(16))*dy(16)));
dy(7)=((k(3)*y(10)*y(1))-(k(4)*y(7)*y(4))+(k(5)*y(5)*y(3))+(k(7)*y(5)*y(6))-(k(8)*y(6))-(k(9)*y(7)*y(2))+ ...
(k(10)*y(10)*y(3))-(k(11)*y(7)*y(9))+(k(12)*y(10)*y(6))-(k(13)*y(7)*(y(10)-y(9)))+(0.5*k(14)*y(6)*(y(11)-y(10)))-(k(15)*y(7)*((y(13)-y(12)))+(0.5*k(16)*y(6)*(y(14)-y(13)))-(0.5*k(20)*(y(8)-y(7)))+(rmon-(y(7)/y(16))*dy(16))));
%mu3
dy(8)=(k(3)*y(11)*y(1))-(k(4)*y(8)*y(4))+(k(5)*y(5)*y(3))+(k(7)*y(5)*((2*y(7))+y(6)))+(k(8)*(y(6)-(2*y(7))))- ...
(k(9)*y(8)*y(2))+(k(10)*y(11)*y(3))-(k(11)*y(8)*y(9))+(k(12)*y(11)*y(6))+(0.33333*k(13)*y(6)*(y(7)-lam3))+ ...
(k(14)*y(7)*(y(8)-y(7)))-(k(15)*y(8)*(y(10)-y(9)))+(0.16667*k(16)*y(6)*((2*mu3))-((3*y(11))+y(10)))-((k(23)*y(8)*(y(13)-y(12))))+(0.166*k(24)*y(6)*((2*gam3)-(3*y(14))+y(13)))-(0.1666*k(20)*((4*lam3)-(3*y(8))-y(7))+(rmon-(y(8)/y(16))*dy(16)));
dy(9)=(-1*k(3)*y(9)*y(1))+(k(4)*y(6)*y(4))+(k(9)*y(6)*y(2))-(k(10)*y(9)*y(3)+(rmon-(y(9)/y(16))*dy(16)));
dy(10)=(-1*k(3)*y(10)*y(1))+(k(4)*y(7)*y(4))+(k(9)*y(7)*y(2))-(k(10)*y(10)*y(3))+(k(11)*y(7)*y(9))-(k(12)*y(10)*y(6))...
+(k(15)*y(7)*(y(10)-y(9)))-(0.5*k(16)*y(6)*(y(11)-y(10)))-(0.5*k(20)*(y(11)-y(10)))+(rmon-(y(10)/y(16))*dy(16));
dy(11)=(-1*k(3)*y(11)*y(1))+(k(4)*y(8)*y(4))+(k(9)*y(8)*y(2))-(k(10)*y(11)*y(3))+(k(11)*y(8)*y(9))-(k(12)*y(11)*y(6))+ ...
(k(15)*y(8)*(y(10)-y(9)))+(k(16)*y(7)*(y(11)-y(10)))+(0.16667*k(16)*y(6)*((-4*mu3)+(3*y(11))+y(10)))-(0.166*k(20)*((4*mu3)-(3*y(11))-y(10))+(rmon-(y(11)/y(16))*dy(16)));
dy(12)=((k(20)*(y(7)-y(6)))+(k(21)*(y(10)-y(9)))+(rmon-((y(12)/y(16))*dy(16))));
dy(13)=(k(13)*y(7)*(y(13)-y(12))-(0.5*k(14)*y(6)*(y(14)-y(13)))-(k(20)*(y(14)-y(13)))...
+(0.5*k(21)*(y(8)-y(7)))+(0.5*k(22)*(y(11)-y(10)))+(rmon-(y(13)/y(16))*dy(16)));
dy(14)=((k(13)*y(8)*(y(13)-y(12))+(k(14)*y(7)*(y(14)-y(13)))+(0.16666*k(15)*y(6)*((-4*gam3+(3*y(14)+y(13)))))-(0.333*k(20)*((4*gam3)-(3*y(14))-y(13))))+(0.1666*k(21)*((2*lam3)-(3*y(8))+y(7)))...
+(0.1666*k(22)*((2*mu3)-(3*y(11))+y(10)))+(rmon-(y(14)/y(16))*dy(16)));
dy(15)=(((((-Hra)*ra)+((-Hrp)*rp))*y(16)))/(mtot*cp);
dy(16)=(((a-(b*diff_mon*dy(15))+c+d+e)*f)-(g*diff_pla*dy(15)));
end
function P=deadkineticestimation(T)
P=zeros(24,1);
ka0=8.11*10^(9);
kp0=8.04*10^(11);
ks0=6.33*10^(7);
kx0=62.1;
kd0=7.29*10^(5);
Ea=29500;
Ep=77700;
Es=9700;
Ex=14500;
Ed=96000;
R=8.314;
Mmon=144.14;
Hra=50100;%%j/mol
Hrp=-23300;
sra=99.3; %j/molK
srp=-22;
P(17)=((1145)/(1+((7.391*10^(-4))*(T-150))));
M0=P(17)/Mmon;
keqa=exp(-1*((Hra-((T+273)*sra)))/(R*(T+273)));
keqp=exp(-1*(Hrp-((T+273)*srp))/(R*(T+273)));
P(1)=ka0*exp((-Ea)/(R*(T+273)));
P(2)=P(1)/keqa;
P(3)=ka0*exp((-Ea)/(R*(T+273)));
P(4)=P(1)/keqa;
P(5)=kp0*exp((-Ep)/(R*(T+273)));
P(6)=((P(3)*M0)/keqp);
P(7)=kp0*exp((-Ep)/(R*(T+273)));
P(8)=((P(3)*M0)/keqp);
P(9)=ks0*exp((-Es)/(R*(T+273)));
P(10)=ks0*exp((-Es)/(R*(T+273)));
P(11)=ks0*exp((-Es)/(R*(T+273)));
P(12)=ks0*exp((-Es)/(R*(T+273)));
P(13)=kx0*exp((-Ex)/((R*(T+273))));
P(14)=kx0*exp((-Ex)/((R*(T+273))));
P(15)=kx0*exp((-Ex)/((R*(T+273))));
P(16)=kx0*exp((-Ex)/((R*(T+273))));
P(18)=exp((-8262.7/(T+273))+6.4361);
P(19)=((-1.4*10^(-2))*(T+273))+7.41;
P(20)=kd0*exp((-Ed)/(R*(T+273)));
P(21)=kd0*exp((-Ed)/(R*(T+273)));
P(22)=kd0*exp((-Ed)/(R*(T+273)));
P(23)=kx0*exp((-Ex)/((R*(T+273))));
P(24)=kx0*exp((-Ex)/((R*(T+273))));
end
% Function call to solve ode
y0=[1 15 0 0 6000 0 0 0 0 0 0 0 0 0 140 10.11];
ts=[0 1];
[P,X]=ode15s(@adia_real,ts,y0,options);
mw=144*(X(:,8)+X(:,11)+X(:,14))./(X(:,10)+(X(:,7)+X(:,13)));
mn=144*(X(:,7)+X(:,10)+X(:,13))./(X(:,9)+X(:,6)+X(:,12));
conv=1-(X(:,5)/y0(5));
pdi=mw./mn ;
plot(P,X(:,15),'-','LineWidth',2.5,'DisplayName','3771');
The plot needs to be like the one shown here...
Thank you
