清空筛选条件
清空筛选条件

Couple projectile motion equations function

2 次查看（过去 30 天）
bml727
bml727 2020-4-28
编辑： James Tursa 2020-12-14
I am creating a function to solve these coupled equations, where all constants are known:
I am currently using the code below, but I am not getting the correct results. All of the constants are specified in my code, b If anyone could offer any help that would be great!
% Set constants for different projectiles (baseball, football, soccer ball)
global m cd A rho g theta v0 x0 y0
m= [0.145 0.42 0.43]; % mass of projectile [kg]
cd= [0.3 0.5 0.25]; % drag coefficients
A= [0.00426 0.02318 0.038]; % cross sectional area of projectile [m^2]
V= [0.0002194 0.00462115 0.00573547]; % volume of projectile
rho= [m(1)/V(1) m(2)/V(2) m(3)/V(3)]; % density
g= 9.8; % gravity [m/s^2]
theta= 45; % initial launch angle [rad]
x0= [0 0 0]; % initial horizontal position
y0= [1.8288 1.8288 0]; % initial vertical position
v0= [18 18 22]; % initial velocity (m/s)
%set projectile, 1=baseball, 2=football, 3=soccer ball
proj=2;
m= m(proj);
cd= cd(proj);
A= A(proj);
rho= rho(proj);
x0= x0(proj);
y0= y0(proj);
v0= v0(proj);
% call function
t_end= 2*v0*sin(theta)/g; % time until projectile reaches the ground
[t y]= ode45(@my_fun,[0 t_end],[x0 y0 v0 theta]);
% plotting results
figure(1)
plot(t,y(:,1))
grid on
xlabel('time [s]')
ylabel('horizontal position [m]')
title('plot of horizontal location vs time')
figure(2)
plot(t,y(:,2))
grid on
xlabel('time [s]')
ylabel('vertical position [m]')
title('plot of vertical location vs time')
figure(3)
plot(y(:,1),y(:,2))
grid on
xlabel('horizontal position [m]')
ylabel('vertical position [m]')
title('plot of vertical location vs horizontal position')
%% Function
function dy=my_fun(t,y) % y is a matrix containing [x y v theta]
global m cd A rho g theta x0 y0 v0
dy=zeros(4,1);
dy(1)= y(3); %dx/dt
dy(2)= y(3); %dy/dt
dy(3)= -(cd*A*rho)/(2*m)*sqrt((dy(1))^2+(dy(2))^2)*dy(1); % d2x/dt2
dy(4)= -g-(cd*A*rho)/(2*m)*sqrt((dy(1))^2+(dy(2))^2)*dy(2); % d2y/dt2
end

请先登录，再回答此问题。

采纳的回答

Ameer Hamza
Ameer Hamza 2020-4-28
In your my_fun you have a line
dy(2)= y(3); %dy/dt
It should be
dy(2)= y(4); %dy/dt
Everything else is correct according to the equation in your question.
  2 个评论
Nicholas Moose
Nicholas Moose 2020-12-14
It looks like y(4) is a launch angle. Shouldn't it be something like dy(1) = y(3)*cos(y(4)) and dy(2) = y(3)*sin(y(4))?
James Tursa
James Tursa 2020-12-14
编辑：James Tursa 2020-12-14
No. There are four states: x, y, vx, vy. These are defined to be y(1), y(2), y(3), and y(4) in the code. The dx/dt value is always going to be vx, which is simply y(3), and the dy/dt value is always going to be vy, which is simply y(4).
The first dx/dt and dy/dt values may in fact correspond to a launch angle, but that is only for the initial point.
That being said, the initial conditions that are passed to ode45( ) do not involve theta directly as was written. They should only be [x0 y0 vx0 vy0] where all of these values are scalars. So in this case since we are working with launch angle theta it would be:
[x0 y0 v0*cosd(theta) v0*sind(theta)]
Also this
theta= 45; % initial launch angle [rad]
should be this
theta= 45; % initial launch angle [deg]
and downstream in the code you should be using sind(theta) instead of sin(theta).

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by