Simulating Projectile with Matlab

2022 3 18
1 个回答
回答已采纳
更新时间：2022 3 20
5 次查看（30 天）

0 个投票

I have a school project where I have to simulate free fall for a basketball from a height of 20m
I have done all the work but when I tried to solve the problem without air resistance (s = (v0 + at ^ 2) / 2) I get that (t =1.3). witch means that the time for the simulation with air resistance is less than the time without air resistance. I dont make sense so I think I have done something wrong with my program but I cant find it. would appreciate the help
thanks

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

 采纳的回答

James Tursa
James Tursa 2022-3-18
编辑：James Tursa 2022-3-18

1 个投票

You don't show us both methods (with and without drag) so we can't compare them. I presume you simply set k=0 for the case without drag? That being said, consider these lines:
V(i+1)=V(i)+deltat*a;
t=(0:N-1)*deltat; % <-- move this out of the loop. It doesn't belong inside the loop.
d=V(i+1)*t + 0.5*a*t.*t % 8 is the starting height
We don't know what your assignment instructions are, but the velocity update you have is an Euler update. If you are supposed to code up Euler, then I would have expected the position update to look similar, e.g. something like this instead of what you have:
d(i+1)=d(i)+deltat*V(i);
And you would start d(1) at some appropriate initial height prior to loop entry. E.g.,
d=zeros(1,N); % Position
d(1)=d0;
As an aside, note that your drag equation only works if the object is always falling down. If it could be going up (e.g., tossed into the air) then the V(i)^2 would have to be replaced with abs(V(i))*V(i). I.e., the sign of the drag effect can change direction depending on which way the ball is moving.
Finally, you should always label all of your constants with units in physical problems like this. It makes it much easier to debug and spot unit problems. Every line that has an explicit number on it should also have the explicit units listed in a comment off to the side.

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Fares Espiro
Fares Espiro 2022-3-18
编辑：Fares Espiro 2022-3-18

Im kind of new to Matlab thats why i dont really know how to do eulers properly. sorry for that
Ye i totaly forgot to do euler with the distance. However i dont understand why it should be:
d(i+1)=d(i)+deltat*V(i);
and not something like:
a = (k/m)*V(i)^2)*t.*t
d(i+1)=d(i)+deltat*V(i)*t + 0.5*(g-a);
My code looks like this now and im getting way better values.
clear all;
V0=0; % initial speed
m=0.615; % mass in kg
g=9.82; % gravity acceleration kg/m3
rho=1.225; % Air density
A=0.04476; % Object area
cw=0.6; % Numerical drag coefficient
k=0.5*cw*rho*A; % Coefficient
N=800; % Time step
V=zeros(1,N); % Speed
V(1)=V0;
deltat=0.002;
t=(0:N-1)*deltat;
d=zeros(1,N);
d(1)=0; % should this be at 0 if im starting from the ground?
for i=1:N-1
a= (g-(k/m)*V(i)^2);
V(i+1)=V(i)+deltat*a;
d(i+1)=d(i)+deltat*V(i)
end
vterminal=sqrt(g*m/k); % Terminal velocity
plot(t,d);
xlabel('time in sec');
ylabel('Distande in m');
legend ('Euler Method','location','south');
XData=get(get(gca,'children'),'YData');
YData=get(get(gca,'children'),'XData');

I was getting around 1.3s without drag and 0.76 with drag
Now im getting around 1.4 with drag.

James Tursa
James Tursa 2022-3-18
"... i dont understand why ..."
The general form for an Euler update of any variable is:
Variable(i+1) = Variable(i) + delta_time * VariableDerivative(i)
If the variable is velocity, then the variable derivative is acceleration.
If the variable is position, then the variable derivative is velocity.
It is as simple as that.
With your code you were apparently trying to combine this simple Euler method with the code for the solution to position based on its 2nd derivative (i.e., the stuff), and my guess is you may have double booked the effect.
Fares Espiro
Fares Espiro 2022-3-18
编辑：Fares Espiro 2022-3-18
Ohh now i see. You are a real hero thank you so much.
Just one last question i want to know the exact time the ball passes d = 8.38
i used this before
interp1(XData,YData,8.38)
but now its giving me an error saying "The grid vectors must contain unique points." how can i fix this?
Fares Espiro
Fares Espiro 2022-3-18
Oh i just figured it out.
I used
[XData, index] = unique(XData);
interp1(XData,YData(index),8.38)
Thank you so much for your help again :)

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

产品

版本

R2019a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by