ode45 how to stop iteration using 'Events'

8 Ansichten (letzte 30 Tage)
Swil
Swil am 6 Dez. 2016
Kommentiert: Walter Roberson am 6 Dez. 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

Antworten (1)

Walter Roberson
Walter Roberson am 6 Dez. 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 Kommentare
Swil
Swil am 6 Dez. 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 am 6 Dez. 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.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by