Unable to meet integration tolerances without reducing the step size

hello, I tried using ODE45 to solve a system of equations however my graph came out looking wonky as well as a warning of [Warning: Failure at t=3.010481e-15. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.310887e-30) at time t.]

clc 
clear
%Forward Reaction
k1 = 4.42E+13;
k2 = 6.91E+12;
k3 = 9.74E+14;
k4 = 8.12E+12;
k5 = 5.61E+10;
k6 = 2.40E+16;
k7 = 7.01E+12;
k8 = 3.64E+11;
k9 = 1.53E+12;
k10 = 6.73+16;
k11 = 9.63+16;
k12 = 4.83E+14;
k13 = 2.72E+09;
k14 = 5.38E+04;
k15 = 4.91E+04;
k16 = 1.08E+08;

Kc1 = 5.20E+04;
Kc2 = 1.81E+04;
Kc5 = 6.21E+04;
Kc9 = 8.6E+04;

tspan = [0 10]';
IC = [10 0 0 0 0 0 0 0 0 0 0 0 0 0 5 5]';
Pt = 1;
Ct0 = .5;
y = 1;
Ft = 10;

% dFdV = @(V,F) [-k1*(((F(1)/Ft)*Ct0*y)-(((F(2)/Ft) *Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc1)), k1*(((F(1)/Ft)*Ct0*y)-(((F(2)/Ft) *Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc1)),k1*(((F(1)/Ft)*Ct0*y)-(((F(2)/Ft) *Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc1))]';
% [V,F_cal] = ode45(dFdV,tspan,IC);
% plot(V,F_cal)
dFdV = @(V,F) [-k1*(((F(1)/Ft)*Ct0*y)-(((F(2)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc1))+k6*((F(9)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y)-k8*((F(1)/Ft)*Ct0*y)^2-k10*((F(2)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y)-k11*((F(4)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y),...%dCa/dT
    k1*(((F(1)/Ft)*Ct0*y)-(((F(2)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc1))-k3*((F(5)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y)-k4*((F(2)/Ft)*Ct0*y)*((F(7)/Ft)*Ct0*y)-k6*((F(9)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y)+k7*((F(4)/Ft)*Ct0*y)^2-k10*((F(2)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y)-k12*((F(2)/Ft)*Ct0*y)^1.24,...%dCb/dT
    k1*(((F(1)/Ft)*Ct0*y)-(((F(2)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc1))+k4*((F(2)/Ft)*Ct0*y)*((F(7)/Ft)*Ct0*y)+k5*(((F(9)/Ft)*Ct0*y)-(((F(4)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc5))+k9*(((F(10)/Ft)*Ct0*y)-(((F(11)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc9)),...%dCc/dT
    -k2*(((F(4)/Ft)*Ct0*y)-(((F(5)/Ft)*Ct0*y)*((F(6)/Ft)*Ct0*y)/Kc2)+k5*(((F(9)/Ft)*Ct0*y)-(((F(4)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc5)))+k6*((F(9)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y)-k7*((F(4)/Ft)*Ct0*y)^2+k10*((F(2)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y)-k11*((F(4)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y)-k13*((F(4)/Ft)*Ct0*y)^1.34,...%dCd/dT
    k2*(((F(4)/Ft)*Ct0*y)-(((F(5)/Ft)*Ct0*y)*((F(6)/Ft)*Ct0*y)/Kc2))-k3*((F(5)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y),...
    k2*(((F(4)/Ft)*Ct0*y)-(((F(5)/Ft)*Ct0*y)*((F(6)/Ft)*Ct0*y)/Kc2))+k8*((F(1)/Ft)*Ct0*y)^2+k10*((F(2)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y),...
    k3*((F(5)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y)-k4*((F(2)/Ft)*Ct0*y)*((F(7)/Ft)*Ct0*y)-k14*((F(7)/Ft)*Ct0*y)^1.37,...
    k4*((F(2)/Ft)*Ct0*y)*((F(7)/Ft)*Ct0*y),...
    -k5*(((F(9)/Ft)*Ct0*y)-(((F(4)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc5))-k6*((F(9)/Ft)*Ct0*y)*((F(2)/Ft)*Ct0*y)+k8*((F(1)/Ft)*Ct0*y)^2,...
    -k9*(((F(10)/Ft)*Ct0*y)-(((F(11)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc9)),...
    k9*(((F(10)/Ft)*Ct0*y)-(((F(11)/Ft)*Ct0*y)*((F(3)/Ft)*Ct0*y)/Kc9))+k11*((F(4)/Ft)*Ct0*y)*((F(1)/Ft)*Ct0*y),...
    k12*((F(2)/Ft)*Ct0*y)^1.24+k13*((F(4)/Ft) *Ct0*y)^1.34+k14*((F(7)/Ft)*Ct0*y)^1.37,...
    -k15*(F(15)*Pt/Ft)-k16*(F(16)*Pt/Ft)^0.31,...
    k15*(F(15)*Pt/Ft)+k16*(F(16)*Pt/Ft)^0.31,...
    -k15*(F(15)*Pt/Ft),...
    -k16*(F(16)*Pt/Ft)^0.31,]';
[V,F_cal] = ode45(dFdV,tspan,IC);
plot(V,real(F_cal))

Answer 1

[prabhat kumar sharma]

Hi Anh,

This problem often arises from stiff equations or rapidly changing solutions where the solver needs to take very small steps to accurately capture the dynamics.

Given your system of equations, there are a few potential issues and solutions you might consider:

1. Check Rate Constants and Initial Conditions

You have k10 = 6.73+16; and k11 = 9.63+16; which seem to be typos. They should likely be k10 = 6.73E+16; and k11 = 9.63E+16;.

2. Adjust ODE Solver Options

You can adjust the solver's error tolerances to make it less strict, although this might result in a less accurate solution. You can do this by setting the AbsTol and RelTol options through odeset.

options = odeset('RelTol',1e-3,'AbsTol',1e-6);
[V,F_cal] = ode45(dFdV,tspan,IC,options);3. Consider Stiff Solvers

3. Consider Stiff Solvers

If your system is stiff, ode45 might not be the best choice as it is designed for non-stiff problems. MATLAB offers solvers like ode15s for stiff problems, which might handle your equations better.

[V,F_cal] = ode15s(dFdV,tspan,IC);

4. Review the Mathematical Model

Rapid changes or discontinuities in the model can cause numerical difficulties. Ensure that the model equations are correctly formulated and that they represent the physical or chemical processes accurately. Pay special attention to the rate equations and any assumptions made.

You can refer this link to understand different solver. https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html#bu3n5rf-1

If the problem persists, you might need to look deeper into the physical model and the numerical methods being used.

I hope it helps!

[Anh Bui]

Thank You so much

[prabhat kumar sharma]

Hi Anh,

If you find this answer helpful , Please accept the answer so that this information can be useful for other community memeber as well.