Unable to meet integration tolerances without reducing the step size
2 views (last 30 days)
Show older comments
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))
0 Comments
Answers (1)
prabhat kumar sharma
on 6 May 2024
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!
2 Comments
prabhat kumar sharma
on 9 May 2024
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.
See Also
Categories
Find more on Ordinary Differential Equations 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!