How to Model Projectile Motion with Air Resistance?

[Maxo]

Grεετings!

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

Homework Statement

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.

Homework 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$$

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:

[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

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?

 

[NascentOxygen]

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

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.

 

[Maxo]

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

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.

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.

 

[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.

 

[Maxo]

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.

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?

 

[NascentOxygen]

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

 

[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?

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.

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?

 

[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.

So the differential equation I'm supposed to solve is a = -k*v ?

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.

 

[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?

 

[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$.

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

Is that correct?

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.

 

[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?

 

[milesyoung]

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? Now I just need to figure out how to write this in MATLAB to make the plots I want.

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.

 

[Maxo]

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.

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?

 

[milesyoung]

Actually I incorporated initial conditions r'(t)=10 for both x and y ...

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.

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?

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.

 

[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?

 

[milesyoung]

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$$

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}$$

 

[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?

 

[milesyoung]

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$$

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.

and that the initial conditions are
$$x(0)=0 \\ x'(0)=10cos(\theta) \\ y(0)=0 \\ y'(0)=10sin(\theta)$$

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

 

[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)?

 

[milesyoung]

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

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.

 

[Maxo]

Okay here are my solutions, I hope you are happy now miles. :thumbs:

$$\\ x(t)=C_1+C_2 e^{-kt} \\ x'(t)=-k C_2 e^{-kt} \\ x(0)=0 \\ x'(0)=10cos \theta \\ x(0)=C_1 + C_2 = 0 \Leftrightarrow C_1 = -C_2 \\ x'(0) = -k C_2 = 10 cos\theta \Leftrightarrow C_2 = -\frac{10}{k} cos \theta \\ \Leftrightarrow C_1 = \frac{10}{k} cos \theta \\ x(t) = \frac{10}{k}cos \theta(1-e^{-kt})$$

$$\\ y(t)=C_3+C_4 e^{-kt}+gt/k \\ y'(t)=-kC_4 e^{-kt}+g/k \\ y(0)=0 \\ y'(0)=sin\theta \\ y(0)=C_3 + C_4 = 0 \Leftrightarrow C_3 = -C_4 \\ y'(0) = -kC_4 + g/k = 10 sin \theta \Leftrightarrow C_4 = -\frac{10}{k} sin \theta + g/k^2 \\ \Leftrightarrow C_3 = g/k^2-\frac{10}{k} sin \theta \\ y(t) = (g/k^2-\frac{10}{k} sin \theta)(1-e^{-kt})+gt/k$$

But they don't give any good plots, so they must still be wrong.

What is wrong?

 

[milesyoung]

First off, let me just say I think you've done a great job. I'm not even sure you were supposed to solve it analytically so, bravo to you sir! :thumbs:

You only have a small error left. You correctly determine that $C_3 = -C_4$, but you don't adhere to it when you write out your expression for $C_3$. As it stands now, $C_3 = C_4$. I think that's what gives you an incorrect expression for $y(t)$. Your $x(t)$ is correct.

 

[Maxo]

First off, let me just say I think you've done a great job. I'm not even sure you were supposed to solve it analytically so, bravo to you sir! :thumbs:

Thank you sir.