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?