Error: Failure at t=1.776273e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.310588e-17) at time t.

I get this error while running my program in ode23s

clc
clear all
close all
t1 = 10;
x0 = [10 0 0 0 0 0 0 0 0 0 0 0 0]';

% for i = 1:length(t1)
    tspan = (0:0.1:t1);
    [t,y] = ode23s(@(t,y) mysystem(t,y),tspan,x0);
    u = y(:,1);
    v = y(:,2);
    w = y(:,3);
    p = y(:,4);
    q = y(:,5);
    r = y(:,6);
    e0 = y(:,7);
    e1 = y(:,8);
    e2 = y(:,9);
    e3 = y(:,10);
    phi = y(:,11);
    theta = y(:,12);
    si = y(:,13);
    x0 = y(end,:);
% end
figure,plot(t,u,'r','linewidth',1.5),hold on,plot(t,v,'--','linewidth',1.5),hold on,plot(t,w,'b','linewidth',1.5)
legend('u','v','w'),grid on;
figure,plot(t,p,'r','linewidth',1.5),hold on,plot(t,q,'--','linewidth',1.5),hold on,plot(t,r,'b','linewidth',1.5)
legend('p','q','r'),grid on;

%%
function dydt =mysystem(t,y)

m = 86816;      
l = 129.5;      
d = 32;         
b = 16;         
Ix = 12470787;  
Iy = 55472495;  
Iz = 49248564;  
Ixz = 1648250;  
ax = 60.39;     
ay = 0.01;      
az = -10.43;   
k1 = 0.085;     
k2 = 0.865;     
k_dash = 0.61;  
Xudot =  -k1*m;
Yvdot = -k2*m;
Zwdot = Yvdot;
Iy_bar = m*(l^2+d^2)/20;
Mqdot = -k_dash*Iy_bar;
Nrdot = Mqdot;
mx = m-Xudot;
my = m-Yvdot;
mz = m-Zwdot;
my = mz;
Jx = Ix;
Jy = Iy-Mqdot;
Jz = Iz-Nrdot;
Jxz = Ixz;
M = [mx 0 0 0 m*az 0;0 my 0 -m*az 0 m*ax;0 0 mz 0 -m*ax 0;0 -m*az 0 Jx 0 -Jxz;m*az 0 -m*ax 0 Jy 0;0 m*ax 0 -Jxz 0 Jz];
mu = 0;
V = 70870.2;    
% S = V^(2/3);    
S = 1;
row = 1.225;
U0 = 25;        
a1 = 60.23;     
a2 = 69.31;     
lf1 = 109.1;
lf2 = 113.7;
lf3 = 13.8;
lh = 103.6;
lgx = 45.3;
lgz = 17.6;
xcv = 63.47;
Sh = 1689.3;
Sf = 709.7;
Sg = 24.7;
etaf = 0.3789;
etak = 1.0518;
CDh0 = 0.0250;
CDf0 = 0.006;
CDg0 = 0.01;
CDch = 0.50;
CDcf = 1;
CDcg = 1;
Ctf = 0;
doalpha = 5.73;
dodelta = 1.24;
f = (lh-a1)/a2;
I1 = pi*((b^2)/(V^(2/3)))*(1-f^2);
I3 = pi*((b^2)/(3*l*V^(2/3)))*(a1-2*a2*f^3-3*a1*f^2)-xcv*I1/l;
J1 = (b/(2*V^(2/3)))*(pi*a1+a2*(sqrt(1-f^2))*f+2*a2*asin(f));
J2 = (J1/l)*(a1-xcv)+((b*2)/(3*l*V^(2/3)))*(a2^2-a1^2-a2^2*((1-f^2)^(2/3)));
Cx1 = -(CDh0*Sh+CDf0*Sf+CDg0*Sg);
Cx2 = (k2-k1)*etak*I1*Sh;
Cx3 = Ctf*Sf;
Cy1 = -Cx2;
Cy2 = -0.5*doalpha*Sf*etaf;
Cy3 = -(CDch*J1*Sh+CDcf*Sf+CDcg*Sg);
Cy4 = -0.5*dodelta*Sf*etaf;
Cz1 = -Cx2;
Cz2 = Cy2;
Cz3 = -(CDch*J1*Sh+CDcf*Sf);
Cz4 = Cy4;
Cl1 = dodelta*Sf*etaf*lf3;
Cl2 = CDcg*Sg*lgz;
Cl3 = CDcg*Sg*lgz*l^2;
Cl4 = -(CDcg*Sg*lgz+CDcf*Sf*lf3)*d^2;
Cm1 = -(k2-k1)*etak*I3*Sh*l;
Cm2 = -0.5*doalpha*Sf*etaf*lf1;
Cm3 = -(CDch*J2*Sh*l+CDcf*Sf*lf2);
Cm4 = -0.5*dodelta*Sf*etaf*lf1;
Cm5 = -CDcf*Sf*lf2*l^2;
Cn1 = -Cm1;
Cn2 = -Cm2;
Cn3 = CDch*J2*Sh*l+CDcf*Sf*lf2+CDcg*Sg*lgx;
Cn4 = -Cm4;
Cn5 = -(CDcf*Sf*lf2+CDcg*Sg*lgx)*l^2;

g = 9.81;       
B = V*row*g;    
H = 1;
W = B+H*g;      
dx = 3.5;   
dy = 7.87;  
dz = 21;    

%% inputs
Tds = 1;
Tdp = 0;
% Tds = [ones(30,1)*1000;ones(30,1)*0;ones(30,1)*0];
% Tdp = [ones(30,1)*1000;ones(30,1)*0;ones(30,1)*0];
deltart = 0;%[ones(Tend,1)*10*pi/180]*0;
deltarb = deltart;
deltael = 0;%zeros(Tend,1);
deltaer = deltael;

% Vtotal = sqrt(y(1)^2+y(2)^2+y(3)^2);
% alpha = atan(y(3)/y(1));
% beta = asin(y(2)/Vtotal);
% 
% Vtotal = sqrt(y(1)^2+y(2)^2+y(3)^2);
% alpha = atan(y(3)/y(1));
% beta = asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)));

dydt = zeros(13,1);
% for i = 1:90
% fval(1) = -mz*y(3)*y(5)+my*y(6)*y(2)+m*(ax*(y(5)^2+y(6)^2)-az*y(6)*y(4))+0.5*row*U0^2*S*(Cx1*(cos(alpha))^2*(cos(beta))^2+Cx2*(sin(2*alpha)*sin(alpha/2)+sin(2*beta)*sin(beta/2)+Cx3))+(Tds(i)+Tdp(i))*cos(mu)+DCM(3,1)*(W-B);
% fval(2) = -mx*y(1)*y(6)+mz*y(4)*y(3)+m*(-ax*y(4)*y(5)-az*y(6)*y(5))+0.5*row*U0^2*S*(Cy1*cos(beta/2)*sin(2*beta)+Cy2*sin(2*beta)+Cy3*sin(beta)*sin(abs(beta))+Cy4*(deltart+deltarb))+DCM(3,2)*(W-B);
% fval(3) = -my*y(2)*y(4)+mx*y(5)*y(1)+m*(-ax*y(6)*y(4)+az*(y(5)^2+y(4)^2))+0.5*row*U0^2*S*(Cz1*cos(alpha/2)*sin(2*alpha)+Cz2*sin(2*alpha)+Cz3*sin(alpha)*sin(abs(alpha))+Cz4*(deltael+deltaer))-(Tdp(i)+Tds(i))*sin(mu)+DCM(3,3)*(W-B);
% fval(4) = -y(6)*y(5)*(Jz-Jy)+Jxz*y(4)*y(5)+m*az*(y(1)*y(6)-y(4)*y(3))+0.5*row*U0^2*S*l*(Cl1*(deltael-deltaer+deltarb-deltart)+Cl2*(sin(beta)*sin(abs(beta)))+0.5*row*S*(Cl3*y(6)*abs(y(6))+Cl4*y(4)*abs(y(4))))+(Tdp(i)-Tds(i))*sin(mu)*dy-DCM(3,2)*az*W;
% fval(5) = -y(4)*y(6)*(Jx-Jz)+Jxz*(y(6)^2-y(4)^2)+m*(ax*(y(2)*y(4)-y(5)*y(1))-az*(y(3)*y(5)-y(6)*y(2)))+0.5*row*U0^2*S*l*(Cm1*cos(alpha/2)*sin(2*alpha)+Cm2*sin(2*alpha)+Cm3*sin(alpha)*sin(abs(alpha))+Cm4*(deltael+deltaer)+0.5*row*S*Cm5*y(5)*abs(y(5)))+(Tds(i)+Tdp(i))*(dz*cos(mu)-dx*sin(mu))+(DCM(3,1)*az-DCM(3,3)*ax)*W;
% fval(6) = -y(5)*y(4)*(Jy-Jx)-Jxz*y(5)*y(6)+m*(-ax*(y(1)*y(6)-y(4)*y(3)))+0.5*row*U0^2*S*l*(Cn1*cos(beta/2)*sin(2*beta)+Cn2*sin(2*beta)+Cn3*sin(beta)*sin(abs(beta))+Cn4*(deltart+deltarb)+0.5*row*S*Cm5*y(6)*abs(y(6)))+(Tdp(i)-Tds(i))*cos(mu)*dy+DCM(3,2)*ax*W;
% end



dydt(1) = -mz*y(3)*y(5)+my*y(6)*y(2)+m*(ax*(y(5)^2+y(6)^2)-az*y(6)*y(4))+0.5*row*U0^2*S*(Cx1*(cos(atan(y(3)/y(1))))^2*(cos(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))^2+Cx2*(sin(2*(atan(y(3)/y(1))))*sin((atan(y(3)/y(1)))/2)+sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))*sin((asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))/2)+Cx3))+(Tds+Tdp)*cos(mu)+2*(y(8)*y(10)-y(7)*y(9))*(W-B);
dydt(2) = -mx*y(1)*y(6)+mz*y(4)*y(3)+m*(-ax*y(4)*y(5)-az*y(6)*y(5))+0.5*row*U0^2*S*(Cy1*cos((asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))/2)*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cy2*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cy3*sin(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))*sin(abs(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cy4*(deltart+deltarb))+2*(y(7)*y(8)-y(9)*y(10))*(W-B);
dydt(3) = -my*y(2)*y(4)+mx*y(5)*y(1)+m*(-ax*y(6)*y(4)+az*(y(5)^2+y(4)^2))+0.5*row*U0^2*S*(Cz1*cos((atan(y(3)/y(1)))/2)*sin(2*(atan(y(3)/y(1))))+Cz2*sin(2*(atan(y(3)/y(1))))+Cz3*sin(atan(y(3)/y(1)))*sin(abs(atan(y(3)/y(1))))+Cz4*(deltael+deltaer))-(Tdp+Tds)*sin(mu)+((y(7)^2)-(y(8)^2)-(y(9)^2)+(y(10)^2))*(W-B);
dydt(4) = -y(6)*y(5)*(Jz-Jy)+Jxz*y(4)*y(5)+m*az*(y(1)*y(6)-y(4)*y(3))+0.5*row*U0^2*S*l*(Cl1*(deltael-deltaer+deltarb-deltart)+Cl2*(sin(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))*sin(abs(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))))+0.5*row*S*(Cl3*y(6)*abs(y(6))+Cl4*y(4)*abs(y(4))))+(Tdp-Tds)*sin(mu)*dy-2*(y(7)*y(8)-y(9)*y(10))*az*W;
dydt(5) = -y(4)*y(6)*(Jx-Jz)+Jxz*(y(6)^2-y(4)^2)+m*(ax*(y(2)*y(4)-y(5)*y(1))-az*(y(3)*y(5)-y(6)*y(2)))+0.5*row*U0^2*S*l*(Cm1*cos((atan(y(3)/y(1)))/2)*sin(2*(atan(y(3)/y(1))))+Cm2*sin(2*(atan(y(3)/y(1))))+Cm3*sin(atan(y(3)/y(1)))*sin(abs(atan(y(3)/y(1))))+Cm4*(deltael+deltaer)+0.5*row*S*Cm5*y(5)*abs(y(5)))+(Tds+Tdp)*(dz*cos(mu)-dx*sin(mu))+2*(y(8)*y(10)-y(7)*y(9))*az-((y(7)^2)-(y(8)^2)-(y(9)^2)+(y(10)^2))*ax*W;
dydt(6) = -y(5)*y(4)*(Jy-Jx)-Jxz*y(5)*y(6)+m*(-ax*(y(1)*y(6)-y(4)*y(3)))+0.5*row*U0^2*S*l*(Cn1*cos((asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))/2)*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cn2*sin(2*(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cn3*sin(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2))))*sin(abs(asin(y(2)/(sqrt(y(1)^2+y(2)^2+y(3)^2)))))+Cn4*(deltart+deltarb)+0.5*row*S*Cm5*y(6)*abs(y(6)))+(Tdp-Tds)*cos(mu)*dy+2*(y(7)*y(8)-y(9)*y(10))*ax*W;
dydt(7) = 0.5*(-y(4)*y(8)-y(5)*y(9)-y(6)*y(10));
dydt(8) = 0.5*(y(4)*y(7)-y(5)*y(10)+y(6)*y(9));
dydt(9) = 0.5*(y(4)*y(10)+y(5)*y(7)-y(6)*y(8));
dydt(10) = 0.5*(-y(4)*y(9)+y(5)*y(8)+y(6)*y(7));
dydt(11) = y(4)+y(5)*tan(y(12))*sin(y(11))+y(6)*tan(y(12))*cos(y(12));
dydt(12) = y(5)*cos(y(11))-y(6)*sin(y(11));
dydt(13) = y(5)*sin(y(11))*sec(y(12))+y(6)*cos(y(11))*sec(y(12));
% fval(7) = 0.5*(-y(4)*y(8)-y(5)*y(9)-y(6)*y(10));
% fval(8) = 0.5*(y(4)*y(7)-y(5)*y(10)+y(6)*y(9));
% fval(9) = 0.5*(y(4)*y(10)+y(5)*y(7)-y(6)*y(8));
% fval(10) = 0.5*(-y(4)*y(9)+y(5)*y(8)+y(6)*y(7));

dydt(1:6) = M\dydt(1:6);
end

Answer 1

[Shadaab Siddiqie]

From my understanding you are getting an error while running above code.

• This is because MATLAB is trying to reduce the time step to a really small amount in order to counter the abrupt change due to the discontinuity in the reference signal.
• For sharp discontinuities,it might not be possible to avoid this warning. However for non discontinous inputs we can set relative and absolute tolerance to a higher number than the default setting.

We can set the tolerances to a higher value for example by using the following commands:

options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode23s(@objectivefunction,tspan,y0,options);