# Projectile with air resistance

1. Mar 29, 2014

### Maxo

Grεετings!

I have some questions about a physics problem. I would appreciate some help with understanding this.

1. The problem statement, all variables and given/known data
The task is to plot some graphs for a projectile at different angles, first when neglecting air resistance and then including a formula for air resistance. The projectile initial speed is 10 m/s and you are supposed to plot the projectile movement from three different initial angles: 20 degrees, 40 degrees and 60 degrees.

2. Relevant equations
We are supposed to use a particular formula for air resistance so it has already been given in the task and it looks like this:
$$\vec{a} = \vec{F}/m = -k\cdot \vec{v} = -k\cdot(v_{x},v_{y})\ with \ k=1$$
3. The attempt at a solution
The plots are supposed to be made with MATLAB. I have already solved the task of plotting the projectile movements when neglecting air resistance, but I need some help when doing it with air resistance.

First of all I wonder if I understand the given formula correctly. I haven't seen a formula for air resistance being written in this way before. When trying to figure it out I would first simplify it, and since k=1 we get simply that $$\vec{a} = -\vec{v}$$ and this would mean that the acceleration in any given point is equal in magnitude to the velocity at that point but opposite in direction. So in the beginning the acceleration would be -10 m/s^2 and in the end it would be +10 m/s^2. Is that a correct interpretation of the formula? I guess this is a possible model to use. But I don't know how to implement it.

I tried doing it in two ways, but none of them worked as I intended.

First try:
Code (Text):
[COLOR="Green"]% Note that variables x and y are vectors of the x and y location of the projectile[/COLOR]
p20 = polyfit(x,y,2);
px = linspace(0,12,100);
fp20 = polyval(p20,px);
% plot(px,fp20,'k'); % this plot is equal to the one I've plotted when excluding air resistance.
v20 = polyder(p20);
fv20 = polyval(v20,px);
fv20x = fv20 * cos(angle(1)); [COLOR="DarkGreen"]% angle(1) being equal to the number of radians corresponding to 20 degrees[/COLOR]
fv20y = fv20 * sin(angle(1));
a20x = -fv20x;
a20y = -fv20y;
x20 = vix .* t + (1/2) * a20x .* t.^2; [COLOR="DarkGreen"]% vix and viy being equal to 10*cos(angle(1)) and 10*sin(angle(1))[/COLOR]
% y20 = viy .* t + (1/2) * (-9.81 - a20y) .* t.^2;
plot(x20,y20,'k');

Second try
Code (Text):
t = linspace(0,2,N);
dh = diff(y);
dt = diff(t);
v = [0 dh./dt];
vx = v * cos(angle(1));
vy = v * cos(angle(1));
a20x = -vx;
a20y = -vy;
x20 = vix .* t + (1/2) * a20x .* t.^2;
y20 = viy .* t + (1/2) * (-9.81 - a20y) .* t.^2;
plot(x20,y20,'k');
The first one is wrong, as it looks (almost) exactly like the plot where air resistance is neglected. The second one I'm not really sure about, since it does give a plot which gives a shorter x-displacement for the 20 and the 40 degree projectile which one would expect, but unfortunately it gives a longer x-displacement for the 60 degree projectile which seems wrong.

Can anyone see what might be wrong?

Last edited: Mar 29, 2014
2. Mar 29, 2014

### Staff: Mentor

Hi maxo. http://img96.imageshack.us/img96/5725/red5e5etimes5e5e45e5e25.gif [Broken]

If there is air resistance, I think the projectile will not return to earth with the same speed as it was launched. You suggested it would, in your dialogue.

Last edited by a moderator: May 6, 2017
3. Mar 30, 2014

### Maxo

That's true. So the acceleration won't be 10 m/s^2 when it hits the ground but just the negative of whatever it's velocity is at that point. I didn't include that assumption in my calculations (MATLAB code) though. It's rather them I need help with now.

If someone could please see what I did wrong in my MATLAB code I would appreciate it. I already saw a misstake myself, that I wrote "vy = v * sin(angle(1));" it's of course supposed to be "vy = v * sin(angle(1));". But that still didn't make it work. Please help.

Last edited by a moderator: May 6, 2017
4. Mar 30, 2014

### milesyoung

You're using a formula that only applies when the acceleration is constant:
http://hyperphysics.phy-astr.gsu.edu/hbase/acons.html

That's fine for the case without air resistance, but you need to solve the governing differential equation from scratch for the case where air resistance makes the acceleration a function of velocity.

This is basically an exercise in calculus unless your assignment just calls for a numerical solution - in that case, you could just have MATLAB supply you with one.

5. Mar 30, 2014

### Maxo

I see. You mean the formula "x = vix .* t + (1/2) * ax .* t.^2"? So that only works for no air resistance? Ok.

I'd be happy to do some calculus to solve this. So the differential equation I'm supposed to solve is a = -k*v ? Which is the same as x'' + k*x' = 0 and y'' + k*y'=0. Right? And then, when I have those solutions, how do I use them?

Last edited: Mar 30, 2014
6. Mar 30, 2014

### Staff: Mentor

Are you sure that you are recomputing the new angle of the velocity after each point? So that you know what acceleration to apply over the next interval? You have that angle(1) term where I was expecting to see a newly calculated angle.

Disclaimer: I don't know Matlab

7. Mar 30, 2014

### Maxo

I did some calculus and I got the solution to the differential equation s(t) = 10(1-e^-t)
(s being displacement). Is that correct?

No I didn't. The angle I uses is constant, but you are right it shouldn't be that. How can I write to retrieve the angle at each given point?

8. Mar 30, 2014

### milesyoung

It looks a bit like, in your MATLAB code, that you want to assume constant acceleration for small increments in time, a la a numerical integration loop, but you're missing the whole 'recalculate stuff in a loop' part.

I think it should be:
$$\mathbf{a} = \mathbf{g} - k \mathbf{v} \Leftrightarrow \frac{\mathrm{d}^2 \mathbf{r}}{{\mathrm{d}t}^2} = \mathbf{g} - k \frac{\mathrm{d} \mathbf{r}}{\mathrm{d}t}$$
where $\mathbf{r}$ is the position vector of the projectile.

Solving it analytically is fairly simple, but you could also just try using, for instance, the 'ode45' command in MATLAB to solve it numerically.

9. Mar 30, 2014

### Maxo

Ok, I solved it and I got the solution:

r(t) = 10-g+(g-10)e^(-k*t)+(g/k)*t

Is that correct?

Now how can I use this solution to solve for the rest of the problem? I need to find a way to get the x and y components of r(t). How can I get that?

10. Mar 31, 2014

### milesyoung

Start, for instance, by decomposing the diff. eq. I wrote into horizontal and vertical components:
\begin{align} x'' &= -k x' &(1) \\ y'' &= g -k y' &(2) \end{align}
You can convert (2) into a first order diff. eq. using the substitution $u = y'$ and solve for $u$ using whatever method you're comfortable with (I used an integrating factor). Solving (2) gives you the solution for (1) as well with $g = 0$.

Some of it is, assuming this is for (2). I can't help further if you don't show step-by-step what you're doing.

11. Mar 31, 2014

### Maxo

Thanks miles, I think I understood the most important things now. I have these solutions:

$$x(t) = 10(1-e^{-kt}) \\ y(t) = 10-g+(g-10)e^{-kt} +\frac{gt}{k}$$
and since k=1 the actual solutions are:
$$x(t) = 10(1-e^{-t}) \\ y(t) = 10-g+(g-10)e^{-t} + gt$$

Right?

But when I plot these two in MATLAB, for which initial angle does the plot show?

Last edited: Mar 31, 2014
12. Mar 31, 2014

### milesyoung

Your solutions still need a bit of work to be correct, but I can't tell you where exactly since you haven't shown their derivations.

It looks like you've tried to incorporate your initial conditions, e.g. $y(0) = 0, \, y'(0) = 10 \sin(20 \cdot \frac{\pi}{180})$, but that still leaves you with some problems.

13. Mar 31, 2014

### Maxo

Actually I incorporated initial conditions r'(t)=10 for both x and y, so I see I made a misstake there. But wouldn't the initial conditions actually be what you wrote, i.e. $$x'(0) = 10 \cos(θ \cdot \frac{\pi}{180}) \\ y'(0) = 10 \sin(θ \cdot \frac{\pi}{180})$$
where θ = 20, 40 or 60 degrees, depending on which case it is?

14. Mar 31, 2014

### milesyoung

You should be careful with what you're doing here, since $\mathbf{r}$ is a vector. I take it you want to specify the magnitude of the velocity vector $\mathbf{r'}$ as an initial condition. That approach will probably only complicate things for you.

I suggest you treat it component-wise, like I showed with (1)-(2), and leave aside the vector part from that point on.

Sure, but have a look at, for instance, your solution $x(t)$:
$$x(t) = 10(1 - e^{-k t}) \Rightarrow x'(t) = 10 k \, e^{-k t} \Rightarrow x'(0) = 10 k$$
So you see, something must be off.

I didn't mean to imply that your x(t) and y(t) aren't solutions to (1) and (2), respectively, but they aren't the solutions (droids) you're looking for in terms of your initial value problem.

I need to see your work to be of any more help to you.

Edit:
I think you're solving your initial value problem correctly, you just aren't using the correct initial conditions when doing it.

Last edited: Mar 31, 2014
15. Mar 31, 2014

### Maxo

Thanks miles. I appreciate your help (and Star Wars reference )

I tried again, and now I got these solutions

$$x(t) = 10k cos \theta (1-e^{kt}) \\ y(t)=(\frac{10sin \theta }{k}-\frac{g}{k})(1-e^{-kt})+\frac{g}{k}t$$

And these seem to give realistic plots for air resistance when k=1. However, according to the assignment, when I you let k → 0 you should get the same plots as the plots for no air resistance. But when I try I don't get those unfortunately. In fact, the lower the value is for k, the smaller the projectile plots become. But they should become the same plots as the ones without air resistance. But they don't. Now, what is wrong?

16. Mar 31, 2014

### milesyoung

Something is still off, but you're getting closer. You can easily check them yourself, e.g.:
\begin{align} x(t) &= 10 k \cos(\theta) (1 - e^{kt}) \Rightarrow x'(t) = -10 k^2 \cos(\theta) e^{kt} \Rightarrow x'(0) = -10 k^2 \cos(\theta)\\ y(t) &= \left(\frac{10 \sin(\theta)}{k} - \frac{g}{k} \right)(1 - e^{-kt}) + \frac{g}{k} t \Rightarrow y'(t) = \left(\frac{10 \sin(\theta)}{k} - \frac{g}{k} \right)(k e^{-kt}) + \frac{g}{k} \Rightarrow y'(0) = 10 \sin(\theta) - g + \frac{g}{k} \end{align}
You should have "gotten your initial conditions back", i.e.:
\begin{align} x'(0) &= 10 \cos(\theta)\\ y'(0) &= 10 \sin(\theta) \end{align}

17. Mar 31, 2014

### Maxo

Do you agree that the general solutions are the following?

$$x(t)=C_{1}+C_{2}e^{-kt} \\ y(t)=C_{1}+C_{2}e^{-kt} + \frac{g}{k}t$$
and that the initial conditions are
$$x(0)=0 \\ x'(0)=10cos(\theta) \\ y(0)=0 \\ y'(0)=10sin(\theta)$$
Is that at least correct?

Last edited: Mar 31, 2014
18. Mar 31, 2014

### milesyoung

Should be:
$$x(t)=C_{1}+C_{2}e^{-kt}\\ y(t)=C_{3}+C_{4}e^{-kt} + \frac{g}{k}t$$
but yes!

Now you just need to figure out what your constants of integration need to be to solve your initial value problem, e.g:
$$x(t) = C_{1}+C_{2}e^{-kt} \Rightarrow x'(t) = -C_2 k e^{kt} \Rightarrow x'(0) = -C_2 k$$
determines $C_2$ and so on.

Looks good. Initial position $(x_0,y_0) = (0,0)$.

19. Mar 31, 2014

### Maxo

I'm starting to get tired and confused, also because Wolfram Alpha gives another answer for y(t):

http://www.wolframalpha.com/input/?i=y''(t)=g-k*y'(t)

As you can see they also divided the "C4"-term with k.

How can that be correct if also what I wrote is correct?

How do you find the particular solution for y(t)?

20. Mar 31, 2014

### milesyoung

Because $c_1$ and $k$ are both constants. Does it really matter if we choose to let $C_4 = \frac{c_1}{k}$? It's the same solution.