ode45 how to stop iteration using 'Events'

11 vues (au cours des 30 derniers jours)
Swil
Swil le 6 Déc 2016
Commenté : Walter Roberson le 6 Déc 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

Réponses (1)

Walter Roberson
Walter Roberson le 6 Déc 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 commentaires
Swil
Swil le 6 Déc 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 le 6 Déc 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.

Connectez-vous pour commenter.

Community Treasure Hunt

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

Start Hunting!

Translated by