Numerical approximation of projectile motion with air resistance

For a school project, I need to estimate the maximum distance of a projectile without neglecting air resistance. I'm following the procedure outlined here:

http://wps.aw.com/wps/media/objects/877/898586/topics/topic01.pdf

On the second page it shows a nice, step by step process to find a numerical approximation. I am trying to reproduce the trajectory of the baseball that is shown on the last page in order to verify my model. However, the plot shows that the baseball will travel over 100 meters, while my model shows that it will travel about 80. Here is my code:

clear 
clc
close all

%Constants for Baseball
  m=.145 %kg
  A=pi*(.0366*2)^2/4 %m^2
  C=.5 %Drag Coefficient of a sphere
  rho= 1.2 %kg/m^3 (density of air)
  D=rho*C*A/2
  g=9.81 %m/s^2 (acceleration due to gravity)

    %Initial Conditions
    delta_t= .001 %s
    x(1)=0
    y(1)=0
    v=50 %m/s
    theta=35 %deg
    vx=v*cosd(theta)
    vy=v*sind(theta)
    t(1)=0 

    %Start Loop
    i=1
    while min(y)> -.001
        ax=-(D/m)*v*vx;
        ay=-g-(D/m)*v*vy;
        vx=vx+ax*delta_t;
        vy=vy+ay*delta_t;

        x(i+1)=x(i)+vx*delta_t+.5*ax*delta_t^2;
        y(i+1)=y(i)+vy*delta_t+.5*ay*delta_t^2;
        t(i+1)=t(i)+delta_t;

        i=i+1;
    end

    plot(x,y)
    xlabel('x distance (m)')
    ylabel('y distance (m)')
    title('Projectile Path')

If anyone has a few minutes to take a look at this, I'd really appreciate it. I can't seem to find the problem. Thank you!

[Walter Roberson]

Well-presented question! You tell us what you are trying to do, you give a reference, you show us your code, you tell us the difference between what you expect and what you get. Thank you!

[Sam]

There is an if statement that has air resistance act in one direction as the ball is travelling upward and act in the other way as the ball travels downward. This is because air is always acting against the motion and the motion changes as the ball travels.

Answer 1

[Teja Muppirala]

You are forgetting to update v inside your loop.

v = sqrt(vx^2 + vy^2);

Right now it is sitting at v = 50 the whole time.

[Dan]

Thanks! that was it

Answer 2

[Stiles]

Also

A=pi*(.0366*2)^2/4 %m^2

should be: A=pi*(.0366)^2 %m^2

as it is the projected surface area

[Stiles]

Updated code showing with and without drag:

clear 
clc
close all
%Constants for Baseball
  m=0.145 %kg
  A=pi*(.0366)^2 %m^2
  C=.5 %Drag Coefficient of a sphere
  rho= 1.2 %kg/m^3 (density of air)
  D=rho*C*A/2
  g=9.81 %m/s^2 (acceleration due to gravity)
    %Initial Conditions
    delta_t= .001 %s
    x(1)=0
    y(1)=0
    x2(1)=0
    y2(1)=0
    v=50 %m/s
    theta=35 %deg
    vx=v*cosd(theta)
    vy=v*sind(theta)
    vx2=v*cosd(theta)
    vy2=v*sind(theta)
    t(1)=0 
    t2(1)=0 
    %Start Loop
    i=1
    while min(y2)> -.001
        ay2=-g;
        ax2=0.0;
        vx2=vx2;
        vy2=vy2+ay2*delta_t;
        x2(i+1)=x2(i)+vx2*delta_t+.5*ax2*delta_t^2;
        y2(i+1)=y2(i)+vy2*delta_t+.5*ay2*delta_t^2;
        t2(i+1)=t2(i)+delta_t;
        i=i+1;
    end
    i=1
    while min(y)> -.001
        ax=-(D/m)*(vx*vx+vy*vy)^0.5*vx;
        ay=-g-(D/m)*(vx*vx+vy*vy)^0.5*vy;
        vx=vx+ax*delta_t;
        vy=vy+ay*delta_t;
        x(i+1)=x(i)+vx*delta_t+.5*ax*delta_t^2;
        y(i+1)=y(i)+vy*delta_t+.5*ay*delta_t^2;
        t(i+1)=t(i)+delta_t;
        i=i+1;
    end

      plot(x,y,'-',x2,y2,'-')
      legend('Linear Drag', 'No Drag')
      legend('show')
      xlabel('x distance (m)')
      ylabel('y distance (m)')
      title('Projectile Path')

[Nikolaj Maack Bielefeld]

Hello mr.

Hats off to you for your work, it inspired me a lot!

However, I made a model that compare your output with the analytical solutions to the differential equations and I found an error in the way you calculate ax and ay. It is wrong to say v*vx or v*vy, it does not make a lot of sense, and I do not understand why you calculate it this way?? Instead I found that when just writing vx^2 and vy^2 the plots come out almost perfectly similar to the analytically calculated. :-) Below is a sample that shows the accuracy of your numerical approximation (after I corrected what I mentioned above) :-)

[darova]

Cant understand why you consider acceleration twice?

vx=vx+ax*delta_t;
x(i+1)=x(i)+vx*delta_t+.5*ax*delta_t^2;

[Nikolaj Maack Bielefeld]

@darova

He calculates each step with constant acceleration. That is only an approximation, but it is a very good one for very small delta_t.

The way I understand his loop:

Line 1&2 (with my own correction from above):

ax3=-(D/m)*vx3^2;
ay3=-g-(D/m)*vy3^2;

First he uses the velocity and the air resistance (a function of velocity, google NASA hehe) to calculate the acceleration component in x- and y-direction. In each direction this is (-D*v^2)/m where v is the x- or y-component of the velocity and minus because it is opposite the movement. In the y-direction gravity is added too. So now we have acceleration. As far as I see you must do this step first in order to begin with how much the projectile is slowed down in the first iteration.

Line 2&3:

vx3=vx3+ax3*delta_t;
vy3=vy3+ay3*delta_t;

These formulas are simply formulas valid for constant acceleration. They calculate the new velocity from the initial velocity and the acceleration from above. What this actually means is that the projectile will not have the initial velocity at any point in this calculation. That's part of the explanation why these iterations always end up a little shorter than the analytical solution, as I've demonstrated above. (The other part is that the iterations always experience more negative acceleration because they are not continuous).

Line 4&5:

x3(i+1)=x3(i)+vx3*delta_t+.5*ax3*delta_t^2;
y3(i+1)=y3(i)+vy3*delta_t+.5*ay3*delta_t^2;

This formula (used for each component), is also simply a standard formula valid for constant acceleration. It is used to calculate the new position of the projectile. As stated above this formula actually never uses the initial velocity. It always calculates the new position from the velocity and acceleration, but the velocity and acceleration are both calculated from the velocity of the previous iteration. But that is actually the best one can do, since a calculation like this can never be continuous.

Now, did this answer your question?