Problem with using ode45 event option

I am trying to simulate a rigid body motion and want to terminate the solver when a particular condition i.e. abs(f_max)=abs(f_req) is satisfied. My main code is:

jw=7e-7;
m=0.0045;
lg=0.015;
g=9.81;
k=0.012;
lm=0.01;
h=9.22e-3;
u=0.3;
M=0.003;
tspan = 0:1/2000:0.15;
theta=pi*30/180;
omega=0;
alpha=(m*g*lg*cos(theta)-k*theta)/(jw+m*lg^2);
Y0=[theta;omega;alpha;0]
options = odeset('RelTol',1e-2,'AbsTol',1e-2,'Events',@(t,y) fric_event(t,y,m,M,g,lg,lm,h,k,u));
[T,Y,te,qe,ie] = ode45(@(t,y) type1(t,y,jw,m,lg,g,k),tspan,Y0,options);
N=(k*Y(:,1)+m*(lm+lg*cos(Y(:,1))).*(g-(Y(:,3))*lg)+m*lg*lm*(Y(:,2).*Y(:,2)).*sin(Y(:,1)))/h;
f_max=u*N;
f_req=((m+M)*g-m*lg*(Y(:,3).*cos(Y(:,1))-Y(:,1).*Y(:,2).*Y(:,2)))/2;
figure
plot(T,abs(f_req),T,abs(f_max),T,0)
legend('abs(f_r)','abs(f_m)')

type 1 function is:

function dy=type1(t,y,jw,m,lg,g,k)
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=(m*g*lg*cos(y(1))-k*y(1))/(jw+m*lg^2);
dy(3)=-(m*g*lg*sin(y(1)).*y(2)+k*y(2))/(jw+m*lg^2);
dy(4)=0;
end

and my events function is

function [value,isterminal,direction] = fric_event(t,y,m,M,g,lg,lm,h,k,u)
  f_req=((m+M)*g-m*lg*(y(3).*cos(y(1))-y(1).*y(2).*y(2)))/2;
  N=(k*y(1)+m*(lm+lg*cos(y(1))).*(g-(y(3))*lg)+m*lg*lm*(y(2).*y(2)).*sin(y(1)))/h;
  f_max=u*N;
  value = abs(f_req)-abs(f_max); % Detect zero of event functions
  isterminal = 1; % Stop the integration
  direction = 0; % Direction
end

After running my code I am getting the following graph of abs(f_max) and abs(f_req) vs time

From the figure it is evident that the condition is satisfied several times but event function is not triggered even once and hence the solution is not terminated. I would like to know where I am going wrong.

Any suggestions would be appreciated. Thank you!

Answer 1

[Mischa Kim]

Sourav, you need to improve the accuracy of the integration. E.g. this seems to work:

 options = odeset('RelTol',1e-7,'AbsTol',1e-12,'Events',@(t,y) fric_event(t,y,m,M,g,lg,lm,h,k,u));