1
$\begingroup$

For a ball of mass $m$, radius $r$, falling freely through air of density $p$, the force of air resistance is proportional to the density of the air, the square of the ball's speed, $v$, and the area of a projection of the ball onto a plane perpendicular to the direction of its motion. Thus, in one dimension, the equation of motion is:

$mx'' = mg - 0.5CApv^2$ , where g is the acceleration due to gravity

If the motion is considered in two dimensions instead, then the following equations describe it:

$mx'' = -0.5CApv_xv$

$my'' = mg - 0.5CApv_yv$

With v = $(v_x,v_y)$ and $v = √(v_x^2 + v_y^2)$

Now I want to solve the two second order differential equations above using the fourth-order Runge-Kutta method in MATLAB. I used the following bolck of code:

SOLUTION IN MATLAB

%NOTES:

%Fouth order Runge-Kutta to solve two systems of second order diff
%equations:
%  x'' = -(1/2CAp/m)v_xV = -Bv_xV ; y'' = g - Bv_yV

%For x: use substitutions: x_1 =x ; x_2 = x_1' = v_x
%for y: use : y_1 = y ; y_2 = y_1' = v_y

%Now we have these two pairs of 1st order equations:

%1) x_1' = x_2 = f(x_1,x_2,y_1,y_2);
%2) x_2' = -Bx_2sqrt(x_2^2 + y_2^2) = g(x_1,x_2,y_1,y_2);
%3) y_1' = y_2 = h(x_1,x_2,y_1,y_2);    
%4) y_2' = g - By_2sqrt(x_2^2 + y_2^2) = q(x_1,x_2,y_1,y_2);

g = 9.81;   %acceleration due to gravity
C = 0.5;    %a constant C, set to 0.5
p = 1.225;  %the density of air (in SI units)
m = 2.47e-03;  %the mass of the ball (SI)
r = 20e-03;    %the radius of the ball (SI)
A = pi*(r^2);  %the projection of the ball onto a normal plane, which is a circle of radius r
B = (C*A*p)/(2*m);  %a combination of all the constants
t_stop = 10;    %the end time(how long I want to integrate the motion)
N = 500;        %number of iterations
dt = t_stop/N;     %time step

x_1(1) = 0;         %initializing the x_1 - arrray
x_2(1) = 7*cos(pi/4);   %initializing the x_2 - arrray
y_1(1) = 0;     %initializing the y_1 - arrray
y_2(1) = 7*sin(pi/4);   %initializing the y_2 - arrray
t(1) = 0;       %the initial time

for n = 1:N;    
    t(n+1) = n*dt;

    k1 = dt*x_2(n);
    m1 = dt*-B*x_2(n)*sqrt(x_2(n)^2 + y_2(n)^2);
    L1 = dt*(y_2(n));
    Z1 = dt*(9.81 - B*y_2(n)*sqrt(x_2(n)^2 + y_2(n)^2));

    k2 = dt*(x_2(n)+m1/2);
    m2 = dt*-B*(x_2(n)+ m1/2)*sqrt((x_2(n)+ m1/2)^2 + (y_2(n)+Z1/2)^2);
    L2 = dt*(y_2(n)+ Z1/2);
    Z2 = dt*(9.81 - B*(y_2(n)+Z1/2)*sqrt((x_2(n)+m1/2)^2 + (y_2(n)+Z1/2)^2));

    k3 = dt*(x_2(n)+m2/2);
    m3 = dt* -B*(x_2(n)+ m2/2)*sqrt((x_2(n)+ m2/2)^2 + (y_2(n)+Z2/2)^2);
    L3 = dt*(y_2(n)+ Z2/2);
    Z3 = dt*(9.81 - B*(y_2(n)+Z2/2)*sqrt((x_2(n)+m2/2)^2 + (y_2(n)+Z2/2)^2));    

    k4 = dt*(x_2(n)+ m3);
    m4 = dt*-B*(x_2(n)+ m3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2);
    L4 = dt*(y_2(n)+ Z3);
    Z4 = dt*(9.81 - B*(y_2(n)+Z3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2));   

    x_1(n+1) = x_1(n)+1/6*(k1 + 2*k2 + 2*k3 + k4);
    x_2(n+1) = x_2(n)+1/6*(m1 + 2*m2 + 2*m3 + m4);
    y_1(n+1) = y_1(n)+1/6*(L1 + 2*L2 + 2*L3 + L4);
    y_2(n+1) = y_2(n)+1/6*(Z1 + 2*Z2 + 2*Z3 + Z4);

end

plot(y_1,x_1)

Problem: When I run the above program, it gives me a rather strange plot that doesn't at all match the motion of a projectile. Can anyone see what the problem is in my code?

$\endgroup$
  • $\begingroup$ Try plotting with: plot(x_1,y_1). As an x-y plot of the motion in a resisting medium this is what I would expect. Better would be subplot(2,2,1);plot(t,x_1);subplot(2,2,2);plot(t,x_2);subplot(2,2,3);plot(t,y_1);subplot(2,2,4);plot(t,y_2). Also note you have selected y as positive downwards, a more normal choice would use y positive upwards which would entail a -ve sign on the mg term. $\endgroup$ – Conrad Turner Sep 29 '15 at 8:36
  • $\begingroup$ The sequence of plots you suggested just gives me the same result, just in four steps. How can you tell that I selected y positive downwards? All I did was to define the initial y and initial x as both zero. $\endgroup$ – 2good4this Sep 29 '15 at 8:44
  • $\begingroup$ Those plots look right for the y positive downwards motion, but to change your code for y positive upwards: use g rather than the hard coded 9.81 in the Z assignments, then set g=-9.81, so gravity acts downwards. $\endgroup$ – Conrad Turner Sep 29 '15 at 8:59
1
$\begingroup$

Your plots look OK except you have defined y as positive downwards by using y"=g + etc. If you want the more usual y positive upwards you need to use y"=-g + etc.

Tidied code with plots of individual components of state against time, and y-positive upwards is:

    %NOTES:

%Fouth order Runge-Kutta to solve two systems of second order diff
%equations:
%  x'' = -(1/2CAp/m)v_xV = -Bv_xV ; y'' = g - Bv_yV

%For x: use substitutions: x_1 =x ; x_2 = x_1' = v_x
%for y: use : y_1 = y ; y_2 = y_1' = v_y

%Now we have these two pairs of 1st order equations:

%1) x_1' = x_2 = f(x_1,x_2,y_1,y_2);
%2) x_2' = -Bx_2sqrt(x_2^2 + y_2^2) = g(x_1,x_2,y_1,y_2);
%3) y_1' = y_2 = h(x_1,x_2,y_1,y_2);    
%4) y_2' = g - By_2sqrt(x_2^2 + y_2^2) = q(x_1,x_2,y_1,y_2);

g = -9.81;   %acceleration due to gravity
C = 0.5;    %a constant C, set to 0.5
p = 1.225;  %the density of air (in SI units)
m = 2.47e-03;  %the mass of the ball (SI)
r = 20e-03;    %the radius of the ball (SI)
A = pi*(r^2);  %the projection of the ball onto a normal plane, which is a circle of radius r
B = (C*A*p)/(2*m);  %a combination of all the constants
t_stop = 10;    %the end time(how long I want to integrate the motion)
N = 500;        %number of iterations
dt = t_stop/N;     %time step

x_1(1) = 0;         %initializing the x_1 - arrray
x_2(1) = 7*cos(pi/4);   %initializing the x_2 - arrray
y_1(1) = 0;     %initializing the y_1 - arrray
y_2(1) = 7*sin(pi/4);   %initializing the y_2 - arrray
t(1) = 0;       %the initial time

for n = 1:N;    
    t(n+1) = n*dt;

    k1 = dt*x_2(n);
    m1 = dt*-B*x_2(n)*sqrt(x_2(n)^2 + y_2(n)^2);
    L1 = dt*(y_2(n));
    Z1 = dt*(g - B*y_2(n)*sqrt(x_2(n)^2 + y_2(n)^2));

    k2 = dt*(x_2(n)+m1/2);
    m2 = dt*-B*(x_2(n)+ m1/2)*sqrt((x_2(n)+ m1/2)^2 + (y_2(n)+Z1/2)^2);
    L2 = dt*(y_2(n)+ Z1/2);
    Z2 = dt*(g - B*(y_2(n)+Z1/2)*sqrt((x_2(n)+m1/2)^2 + (y_2(n)+Z1/2)^2));

    k3 = dt*(x_2(n)+m2/2);
    m3 = dt* -B*(x_2(n)+ m2/2)*sqrt((x_2(n)+ m2/2)^2 + (y_2(n)+Z2/2)^2);
    L3 = dt*(y_2(n)+ Z2/2);
    Z3 = dt*(g - B*(y_2(n)+Z2/2)*sqrt((x_2(n)+m2/2)^2 + (y_2(n)+Z2/2)^2));    

    k4 = dt*(x_2(n)+ m3);
    m4 = dt*-B*(x_2(n)+ m3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2);
    L4 = dt*(y_2(n)+ Z3);
    Z4 = dt*(g - B*(y_2(n)+Z3)*sqrt((x_2(n)+ m3)^2 + (y_2(n)+Z3)^2));   

    x_1(n+1) = x_1(n)+1/6*(k1 + 2*k2 + 2*k3 + k4);
    x_2(n+1) = x_2(n)+1/6*(m1 + 2*m2 + 2*m3 + m4);
    y_1(n+1) = y_1(n)+1/6*(L1 + 2*L2 + 2*L3 + L4);
    y_2(n+1) = y_2(n)+1/6*(Z1 + 2*Z2 + 2*Z3 + Z4);

end

subplot(2,2,1);plot(t,x_1);title('x against time');
subplot(2,2,2);plot(t,x_2);title('v_x against time');

subplot(2,2,3);plot(t,y_1);title('y against time');
subplot(2,2,4);plot(t,y_2);title('v_y against time');

Which produces plots:

enter image description here

from which we see that v_x drops to zero, and v_y approaches a terminal velocity.

$\endgroup$
  • $\begingroup$ Brilliant! Thank you! It works perfect now. $\endgroup$ – 2good4this Sep 29 '15 at 9:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.