problems with integration tolerances when using ode15s

I am working on simulating a dynamic system using a Concurrent Learning controller, but for now all I am doing is plotting the norm of my errors (e and r) and the norm of one of my estimates (utilde). When I run my code I get a yellow warning stating "Warning: Failure at t=1.879973e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-14) at time t." I'm not sure if this is due to my reltol and abstol values or a problem within my code but it is effecting how my plots should look. The three plots should show convergence to zero, but norm utilde(t) plots strangely while norm e and norm r seem to plot normally. My code is posted below. Any help of suggestions are greatly appreciated.

clc;
clear;
close all;
low=0;%stated that all constants mx,mphi,cx,cphi,1,ax,and aphi were greater than zero 
high=1;
Amx=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant mx
Amphi=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant mphi
Acx=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant cx
Acphi=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant cphi
Al=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant l
Aax=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant ax
Aaphi=(high-low).*rand(100,1) + low;%column vector of 100 uniformly distributed decimals between 0 and 1 for constant aphi

%Feedback terms
lowc=0;
highc=5;
ABx=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Bx
ABphi=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Bphi
Aalphax=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant alphax
Aalphaphi=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant alphaphi
AJx=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Jx
AJphi=(highc-lowc).*rand(100,1) + lowc;%column vector of 100 uniformly distributed decimals between 0 and 5 for constant Jphi

%Creating xd and phid equations and their derivatives
lowxd=-1;%left limit of xd set to 1m
highxd=1;%right limit of xd set to 1m
Axdbar=(highxd-lowxd).*rand(100,1)+lowxd;%column vector of 100 uniformly distributed decimals between -1 and 1 for variable xdbar
Ax=(highxd-lowxd).*rand(100,1)+lowxd;%column vector of 100 uniformly distributed decimals between -1 and 1 for variable x 
lowfreq=0;%lowest frequency set to 0 Hz
highfreq=0.5;%highest frequency set to 0.5 Hz
Afxd=(highfreq-lowfreq).*rand(100,1)+lowfreq;%column vector of 100 uniformlay distributed decimals between 0 and 0.5 Hz for variable fxd
lowphi=-1*(pi/2);
highphi=pi/2;
Aphidbar=(highphi-lowphi).*rand(100,1)+lowphi;%column vector of 100 uniformly distributed decimals between -pi/2 and pi/2 for variable phidbar
Aphi=(highphi-lowphi).*rand(100,1)+lowphi;%column vector of 100 uniformly distributed decimals between -pi/2 and pi/2 for variable phi 
Afphid=(highfreq-lowfreq).*rand(100,1)+lowfreq;%column vector of 100 uniformlay distributed decimals between 0 and 0.5 Hz for variable fphid

for ii=1:100 %signifies 100 iterations of each constant and variable
    %Constants
    mx=Amx(ii);
    mphi=Amphi(ii);
    cx=Acx(ii);
    cphi=Acphi(ii);
    l=Al(ii);
    ax=Aax(ii);
    aphi=Aaphi(ii);
    
    %Trajectory Variables
    xdbar=Axdbar(ii);
    x=Ax(ii);
    fxd=Afxd(ii);
    phidbar=Aphidbar(ii);
    phi=Aphi(ii);
    fphid=Afphid(ii);
    g=9.8;%gravitational constant
    
    %Feedback Terms
    Bx=ABx(ii);
    Bphi=ABphi(ii);
    alphax=Aalphax(ii);
    alphaphi=Aalphaphi(ii);
    Jx=AJx(ii);
    Jphi=AJphi(ii);
    
    gammat=eye(6);
    gammaa=eye(2);
    Ktheta=eye(6);
    Ka=eye(2);
    
    tspan= [0 50];%time span for elapsed simulation
    options = odeset('reltol', 1e-3, 'abstol', 1e-5 );%default options for relative and absolute tolerance.  The tolerances are used to limit the local discretization error during integration
    x0=[x 0 0 phi 0 0 0 0 0 0 0 0 0 0 0];
    [t, x] = ode15s(@(t,x) dynamics(t,x,phi,mx,mphi,cx,cphi,l,ax,aphi,Bx,Bphi,alphax,alphaphi,Jx,Jphi,xdbar,fxd,phidbar,fphid,g,gammat,gammaa,Ktheta,Ka), tspan, x0, options);
    
    t=transpose(t);
    x=transpose(x);
    x1=x(1,:);
    xdot=x(2,:);
    xdotdot=x(3,:);
    phi1=x(4,:);
    phidot=x(5,:);
    ux=x(6,:);
    uphi=x(7,:);
    thetahat=x(8:13,:);
    ahat=x(14:15,:);
    
    %Trajectory equations
    xd=xdbar*sin(2*pi*fxd.*t);%xd(t)
    xddot=2*pi*fxd*xdbar*cos(2*pi*fxd.*t);%xd(t)dot
    xddotdot=-1*(2*pi*fxd)^(2)*xdbar*sin(2*pi*fxd.*t);%xd(t)dotdot
    phid=phidbar*sin(2*pi*fphid.*t);%phid(t)
    phiddot=2*pi*fphid*phidbar*cos(2*pi*fphid.*t);%phid(t)dot
    phiddotdot=-1*(2*pi*fphid)^(2)*phidbar*sin(2*pi*fphid.*t);%phid(t)dotdot

    %Estimation Errors
    sizet=size(t,2);
    theta=[mx*ones(1,sizet); mphi*ones(1,sizet); mphi*l*ones(1,sizet); mphi*(l)^(2)*ones(1,sizet); cx*ones(1,sizet); cphi*ones(1,sizet)];
    a=[ax*ones(1,sizet); aphi*ones(1,sizet)];
    thetatilde=theta-thetahat;
    atilde=a-ahat;
    
    %Tracking Errors
    ex=xd-x1;
    exdot=xddot-xdot;
    exdotdot=xddotdot-xdotdot;
    ephi=phid-phi1;
    ephidot=phiddot-phidot;
    ephidotdot=phiddotdot-((-1*cphi.*phidot-mphi*l*g*sin(phi)-mphi*l*cos(phi).*xdotdot+uphi)/(mphi*l^(2)));
    
    %Filtered Tracking Errors
    rx=exdot+alphax.*ex;
    rxdot=exdotdot+alphax.*exdot;
    rphi=ephidot+alphaphi.*ephi;
    rphidot=ephidotdot+alphaphi.*ephidot;
    
    %Design Inputs
    yxthetahat=zeros(1,length(t));
    yphithetahat=zeros(1,length(t));
    Yxthetahat=zeros(1,length(t));
    Yphithetahat=zeros(1,length(t));
    yxathetahatdot=zeros(1,length(t));
    yphiathetahatdot=zeros(1,length(t));
    Yxthetahatdot=zeros(1,length(t));
    Yphithetahatdot=zeros(1,length(t));
    for jj=1:length(t)
    yx=[xdotdot(1,jj) (sin(phi))^(2)*xdotdot(1,jj)-g*cos(phi)*sin(phi) phidot(1,jj)^(2)*sin(phi) 0 xdot(1,jj) 0];
    yphi=[0 0 g*sin(phi)+cos(phi)*xdotdot(1,jj) (-1*cphi*phidot(1,jj)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*xdotdot(1,jj)+uphi(1,jj))/(mphi*l^(2)) 0 phidot(1,jj)];   
    Yx=[xddotdot(1,jj)+alphax*exdot(1,jj) sin(phi)^(2)*xdotdot(1,jj)-g*cos(phi)*sin(phi)+sin(phi)^(2)*alphax*exdot(1,jj)+sin(phi)*cos(phi)*phidot(1,jj)*rx(1,jj) phidot(1,jj)^(2)*sin(phi) 0 xdot(1,jj) 0];
    Yphi=[0 0 g*sin(phi)+cos(phi)*xdotdot(1,jj) phiddotdot(1,jj)+alphaphi*ephidot(1,jj) 0 phidot(1,jj)];
    yxa=[ux(1,jj) 0];
    yphia=[0 uphi(1,jj)];
%     thetahatdot=gammat*Yx.'*rx+gammat*Yphi.'*rphi+gammat*Ktheta*(yx.'*x(6)-yx.'*yx*x(8:13))+gammat*Ktheta*(yphi.'*x(7)-yphi.'*yphi*x(8:13));
    yxthetahat(jj)=yx*thetahat(:,jj);
    yphithetahat(jj)=yphi*thetahat(:,jj);
    Yxthetahat(jj)=Yx*thetahat(:,jj);
    Yphithetahat(jj)=Yphi*thetahat(:,jj);
%     Yxthetahatdot(jj)=Yx*thetahatdot
    end
    udx=Yxthetahat+ex+Bx.*rx;
    udphi=Yphithetahat+ephi+Bphi.*rphi;
    utildex=udx-ux;
    utildephi=udphi-uphi;
    
    
%     udxdot=Yx*thetahatdot+exdot+Bx*rxdot;
%     udphidot=Yphi*thetahatdot+ephidot+Bphi*rphidot;
%     mux=udxdot+x(14)*x(6)+rx+Jx*utildex;
%     muphi=udphidot+x(15)*x(7)+rphi+Jphi*utildephi;
%     ahatdot=gammaa*utildex*x(6)+gammaa*utildephi*x(7)+gammaa*Ka*(yxa.'*(mux-(-1*ax*x(6)+mux))-yxa.'*yxa*x(14:15))+gammaa*Ka*(yphia.'*(muphi-(-1*aphi*x(7)+muphi))-yphia.'*yphia*x(14:15));
%     
    %Tracking Error Plot
    norme=sqrt(ex.^2+ephi.^2);
    
    %Filtered Tracking Error Plot
    normr=sqrt(rx.^2+rphi.^2);
        
    %Estimation Errors
    normutilde=sqrt(utildex.^2+utildephi.^2);
    
    figure(1)
    plot(t,norme)
    xlabel('Time')
    ylabel('norm e(t)')
    hold on
    figure(2)
    plot(t,normr)
    xlabel('Time')
    ylabel('norm r(t)')
    hold on
    figure(3)
    plot(t, normutilde)
    xlabel('Time')
    ylabel('norm utilde(t)')
    hold on
%     figure(4)
%     plot(t,normthetatilde)
%     xlabel('Time')
%     ylabel('norm thetatilde(t)')
%     hold on
%     figure(5)
%     plot(t,normatilde)
%     xlabel('Time')
%     ylabel('norm atilde(t)')
%     hold on

end
hold off
function dstate = dynamics(t,x,phi,mx,mphi,cx,cphi,l,ax,aphi,Bx,Bphi,alphax,alphaphi,Jx,Jphi,xdbar,fxd,phidbar,fphid,g,gammat,gammaa,Ktheta,Ka)
    %Trajectory equations
    xd=xdbar*sin(2*pi*fxd*t);%xd(t)
    xddot=2*pi*fxd*xdbar*cos(2*pi*fxd*t);%xd(t)dot
    xddotdot=-1*(2*pi*fxd)^(2)*xdbar*sin(2*pi*fxd*t);%xd(t)dotdot
    phid=phidbar*sin(2*pi*fphid*t);%phid(t)
    phiddot=2*pi*fphid*phidbar*cos(2*pi*fphid*t);%phid(t)dot
    phiddotdot=-1*(2*pi*fphid)^(2)*phidbar*sin(2*pi*fphid*t);%phid(t)dotdot

    %Estimation Errors
    theta=[mx;mphi;mphi*l;mphi*(l)^(2);cx;cphi];
    a=[ax;aphi];
    thetatilde=theta-x(8:13);
    atilde=a-x(14:15);
    
    %Tracking Errors
    ex=xd-x(1);
    exdot=xddot-x(2);
    exdotdot=xddotdot-x(3);
    ephi=phid-x(4);
    ephidot=phiddot-x(5);
    ephidotdot=phiddotdot-((-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2)));
    
    %Filtered Tracking Errors
    rx=exdot+alphax*ex;
    rxdot=exdotdot+alphax*exdot;
    rphi=ephidot+alphaphi*ephi;
    rphidot=ephidotdot+alphaphi*ephidot;
    
    %Design Inputs
    yx=[x(3) (sin(phi))^(2)*x(3)-g*cos(phi)*sin(phi) x(5)^(2)*sin(phi) 0 x(2) 0];
    yphi=[0 0 g*sin(phi)+cos(phi)*x(3) (-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2)) 0 x(5)];   
    Yx=[xddotdot+alphax*exdot sin(phi)^(2)*x(3)-g*cos(phi)*sin(phi)+sin(phi)^(2)*alphax*exdot+sin(phi)*cos(phi)*x(5)*rx x(5)^(2)*sin(phi) 0 x(2) 0];
    Yphi=[0 0 g*sin(phi)+cos(phi)*x(3) phiddotdot+alphaphi*ephidot 0 x(5)];
    udx=Yx*x(8:13)+ex+Bx*rx;
    udphi=Yphi*x(8:13)+ephi+Bphi*rphi;
    yxa=[x(6) 0];
    yphia=[0 x(7)];
    utildex=udx-x(6);
    utildephi=udphi-x(7);
    thetahatdot=gammat*Yx.'*rx+gammat*Yphi.'*rphi+gammat*Ktheta*(yx.'*x(6)-yx.'*yx*x(8:13))+gammat*Ktheta*(yphi.'*x(7)-yphi.'*yphi*x(8:13));
    udxdot=Yx*thetahatdot+exdot+Bx*rxdot;
    udphidot=Yphi*thetahatdot+ephidot+Bphi*rphidot;
    mux=udxdot+x(14)*x(6)+rx+Jx*utildex;
    muphi=udphidot+x(15)*x(7)+rphi+Jphi*utildephi;
    ahatdot=gammaa*utildex*yxa.'+gammaa*utildephi*yphia.'+gammaa*Ka*(yxa.'*(mux-(-1*ax*x(6)+mux))-yxa.'*yxa*x(14:15))+gammaa*Ka*(yphia.'*(muphi-(-1*aphi*x(7)+muphi))-yphia.'*yphia*x(14:15));

    %Define State Variables
    dstate = zeros(15,1);%initialize system matrix
    dstate(1)=x(2);%xdot
    dstate(2)=x(3);%xdotdot
    dstate(3)=(-1*cx*x(3)+mphi*g*x(5)*(cos(phi)^(2)-sin(phi)^(2))+mphi*l*x(5)*(2*((-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2)))*sin(phi)+x(5)^(2)*cos(phi))+(-1*ax*x(6)+(mux))-2*mphi*cos(phi)*x(5)*x(3))/(mx+mphi*sin(phi)^(2));%xdotdotdot
    dstate(4)=x(5);%phidot
    dstate(5)=(-1*cphi*x(5)-mphi*l*g*sin(phi)-mphi*l*cos(phi)*x(3)+x(7))/(mphi*l^(2));%phidotdot
    dstate(6)=-1*ax*x(6)+(mux);%uxdot
    dstate(7)=-1*aphi*x(7)+(muphi);%uphidot
    dstate(8:13)=thetahatdot;%thetahatdot
    dstate(14:15)=ahatdot;%ahatdot
end

Answer 1

[Sulaymon Eshkabilov]

The answer lies in magnitudes of the calculated (integration) values of the variables that don't allow to have the set tolerances without reducing the time step. Therefore, it is better to loosen the tolerances, i.e. reltol and abstol. As it is now, the tolerances are quite significantly tight.

[Matthew Tortorella]

Thank you Sulaymon, that worked. The program still takes about 2-3 minutes to plot all of the iterations but I don't see any integration tolernace issues any more.

[Sulaymon Eshkabilov]

Most Welcome!