Using fourth order Runge-Kutta to solve a second order ode in MATLAB

Question

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?

Answer 1

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:

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