Projectile motion with friction in MATLAB (ODE45)

[cwrn]

I'm working on a little project where I want to plot the motion of a projectile with air resistance. The air resistance can be assumed to be proportional to the velocity squared.

F_{f}=-Bv^{2}

F_{f,x}=F_{f}\frac{v_{x}}{v}, \ \ F_{f,y}=F_{f}\frac{v_{y}}{v}

where B depends on the height, y above the ground

B(y) = B_{0}e^{-y/y_{0}}

Given

k = \frac{B_{0}}{m}=4\cdot 10^{-5} \ m^{-1}, \ y_{0} = 1\cdot 10^4 \ m, \ v_{0} = 700 \ m/s

I have derived the equations of motion as following:

\vec{r}=\left(v_{0}tcos\theta-\frac{1}{2}kvv_{x}e^{-y/y_{0}}\cdot t^2\right)\hat{i}+\left(v_{0}tcos\theta+\frac{1}{2}\left[\vec{g}-kvv_{y}e^{-y/y_{0}}\right] t^2\right)\hat{j}

\vec{v}=\left(v_{0}cos\theta-kvv_{x}e^{-y/y_{0}}\cdot t\right)\hat{i}+\left(v_{0}sin\theta+\left [\vec{g}-kvv_{y}e^{-y/y_{0}}\right]t\right)\hat{j}

\vec{a}=\left(-kvv_{x}e^{-y/y_{0}}\right)\hat{i}+\left(\vec{g}-kvv_{y}e^{-y/y_{0}}\right)\hat{j}

I'm having trouble defining a function I can use with ode45 since there are several variables depending on each other (assuming my equations are correct). Any tips would be greatly appreciated.

 

[DrClaude]

I haven't checked if the equations are correct, but they are not differential equations!

What you need for ode45 is something like
$$
\begin{align}
\frac{dx}{dt} &= v_x \\
\frac{dy}{dt} &= v_y \\
\frac{dv_x}{dt} &= a_x \\
\frac{dv_y}{dt} &= a_y
\end{align}
$$
where you replace $a_x$ and $a_y$ by the proper expressions involving the forces.

 

[cwrn]

I haven't checked if the equations are correct, but they are not differential equations!

What you need for ode45 is something like
$$
\begin{align}
\frac{dx}{dt} &= v_x \\
\frac{dy}{dt} &= v_y \\
\frac{dv_x}{dt} &= a_x \\
\frac{dv_y}{dt} &= a_y
\end{align}
$$
where you replace $a_x$ and $a_y$ by the proper expressions involving the forces.

Yes, I'm aware of this. Solving this with ode45 should yield the x and y vectors that I want to plot. I'm just not sure how to define the function itself since it depends on B which in turn depends on y.
$$
\begin{align}
B(y) = B_{0}e^{-y/y_{0}}
\end{align}
$$

 

[DrClaude]

If you take the array y to contain $(x,y,v_x,v_y)$, then you would write something like

function dy = projectile(t,y)
B0 = 1;
y0 = 123;

dy = zeros(4,1);    % a column vector
dy(1) = y(3);
dy(2) = y(4);
dy(3) = B0 * exp(-y(2)/y0) * y(3);
dy(4) = B0 * exp(-y(2)/y0) * y(4);

(Note that the code doesn't correspond to your problem.)

 

[SteamKing]

Yes, I'm aware of this. Solving this with ode45 should yield the x and y vectors that I want to plot. I'm just not sure how to define the function itself since it depends on B which in turn depends on y.
$$
\begin{align}
B(y) = B_{0}e^{-y/y_{0}}
\end{align}
$$

You start off your calculations with a known position for the projectile, which means you are going to know an initial value for y. As the calculation progresses, then the values of y will need to be updated depending on your independent variable, which I assume would be time.