Simulating Projectile with Matlab

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

Answer 1

[James Tursa]

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.

[Fares Espiro]

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]

"... 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]

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]

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 :)