ode45 how to stop iteration using 'Events'

11 views (last 30 days)
Swil
Swil on 6 Dec 2016
Commented: Walter Roberson on 6 Dec 2016
I have a function of the trajectory of a smoothball and am trying to get the iteration to stop once the ball lands and am having trouble implementing 'Events'. This is my function could any one help with the events function? Thanks!
function Smoothball(theta)
% given values
m = 0.045; % [kg] mass
D = 0.043; % [m] diameter
rho = 1.23; % [kg/m^3] air density
nu = 1.46e-5; % [m^2/s] air kinematic viscosity
g = 9.807; % [m/s^2] gravity
V0 = 80; % [m/s] initial speed
%calculate cross sectional area
A = pi * (D/2)^2; %cross sectional area
[t,x]=ode45(@flights,[0,6.299],[0,V0*cos(theta),0,V0*sin(theta),0,0]);
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8,'Events',@events);
plot(x(:,1),x(:,3));
title('Smooth Ball Model')
xlabel('X')
ylabel('Y')
grid on
end

Answers (1)

Walter Roberson
Walter Roberson on 6 Dec 2016
function [value, isterminal, direction] = events(t,y)
value = y;
isterminal = 1;
direction = 0;
This could also be written as the more obscure
events = @(t, y) deal(y, 1, 0);
in which case the 'Events' parameter would name events rather than @events
  2 Comments
Swil
Swil on 6 Dec 2016
function x = Smoothball(theta)
% given values
m = 0.045; % [kg] mass
D = 0.043; % [m] diameter
rho = 1.23; % [kg/m^3] air density
nu = 1.46e-5; % [m^2/s] air kinematic viscosity
g = 9.807; % [m/s^2] gravity
V0 = 80; % [m/s] initial speed
%calculate cross sectional area
A = pi * (D/2)^2; %cross sectional area
[t,x]=ode45(@flights,[0,6.299],[0,V0*cos(theta),0,V0*sin(theta),0,0])
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8,'Events',@events);
plot(x(:,1),x(:,3));
title('Smooth Ball Model')
xlabel('X')
ylabel('Y')
grid on
end
function xprime=flights(t,x)
%Parameters
C_d=.25;
r=.002378;
A=.25*pi*(1.75/12)^2;
m=.045;
%D =((1/2)*C_d*r*A*V^2
%Square of the velocity is in the equation, therefore
D =((1/2)*C_d*r*A);
xprime=zeros(2,1);
%X
xprime(1)=x(2);
xprime(2)=-(D/m)*x(2)^2;
%Y
xprime(3)=x(4);
xprime(4)=-32.2-(D/m)*x(4)^2;
%Z
xprime(5)=x(6);
xprime(6)=-(D/m)*x(6)^2;
end
function [value, isterminal, direction] = events(t,x)
value = x;
isterminal = 1;
direction = 0;
end
so this is my code now however the iteration is still going into negative y values
Walter Roberson
Walter Roberson on 6 Dec 2016
In events()
value = x(3);
if it is the y you intend to test.
But your larger problem was that you create your options after you call ode45() and you do not pass the options to ode45. You need to set the options first and pass the options structure to ode45.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!